A single-cell transcriptome atlas of the aging human and macaque retina

ABSTRACT The human retina is a complex neural tissue that detects light and sends visual information to the brain. However, the molecular and cellular processes that underlie aging primate retina remain unclear. Here, we provide a comprehensive transcriptomic atlas based on 119 520 single cells of the foveal and peripheral retina of humans and macaques covering different ages. The molecular features of retinal cells differed between the two species, suggesting distinct regional and species specializations of the human and macaque retinae. In addition, human retinal aging occurred in a region- and cell-type-specific manner. Aging of human retina exhibited a foveal to peripheral gradient. MYO9A− rods and a horizontal cell subtype were greatly reduced in aging retina, indicating their vulnerability to aging. Moreover, we generated a dataset showing the cell-type- and region-specific gene expression associated with 55 types of human retinal disease, which provides a foundation to understanding of the molecular and cellular mechanisms underlying human retinal diseases. Such datasets are valuable to understanding of the molecular characteristics of primate retina, as well as molecular regulation of aging progression and related diseases.

(Continued on next page)

INTRODUCTION
The human retina is a specialized light-sensitive tissue of neurons, glia and nourishing blood vessels [1][2][3]. The retina has long served as a model system for developmental and functional studies of the central nervous system [4,5]. Different cell types in the retina, including rod and cone photoreceptors as well as bipolar (BCs), amacrine (ACs), horizontal (HCs), ganglion (GCs) and glial cells, are packed together into a tightly organized network that converts incoming light into electrochemical signals, which are then relayed to the brain for visual formation [1]. In addition to this complexity, primates, including humans and monkeys, possess a specialized fovea, which is absent in rodent models. The fovea, which is responsible for high-acuity vision, is enriched in cone photoreceptors that directly receive light and are supported by specific morphological Müller glia (MGs) [6]. The molecular differences between the foveal and peripheral retina are important for understanding human visual function. To date, multiple analyses in primates have demonstrated the region-specific retinal transcriptomes in developing and adult foveal and peripheral retina, suggesting distinct transcriptional regulations in the two regions [3,[7][8][9][10][11][12].
As people age, a deterioration in photoreceptor structure occurs in the human retina, particularly in the foveal region. Patients aged ≥60 years display obvious reductions in the number of foveal photoreceptors, and a decline in color discrimination and vision sensitivity [13]. Moreover, they are at detailed molecular and cellular signatures between the different cell types in human and macaque retina remain elusive. In particular, whether various cell types in the human and macaque retina exhibit different molecular changes following aging is largely unknown.
As the primate retina is a heterogeneous tissue containing various cell types, the transcriptional heterogeneity of retinal cells cannot be detected using conventional bulk RNA-sequencing (RNA-seq). In contrast, single cell RNA-sequencing (scRNA-seq) has the power to detect region-and cell-type-specific changes, especially for rare cell types. In addition, scRNA-seq can decipher ligand-receptor crosstalk between cell-cell interactions [18]. To decipher the detailed molecular processes that accompany retinal aging, we mapped >119,000 cell transcriptomes of the human and macaque (Macaca mulatta) retina from young to old stages and explored related gene transcriptional regulation using scRNA-seq and the bulk Assay for Transposase-Accessible Chromatin with high-throughput sequencing (bulk ATAC-seq). Our data revealed that although the human and macaque retina showed high similarity in major cell types, the molecular features were distinguishable between species. Rods showed the greatest interspecies variation. Although rods from humans and macaques could be divided into two subtypes by MYO9A expression, the proportion of subtypes varied between the species. Interestingly, MGs and cones not only exhibited regional transcription differences, they also interacted with each other by ligand-receptor regulation. Furthermore, aging in the human foveal retina occurred earlier than that of the peripheral retina, and MYO9A − rods and one HC subtype showed greater vulnerability to aging. Finally, we generated a dataset showing the cell-type-and region-specific gene expression associated with 55 types of human retinal disease. Overall, our study provides a valuable data source for understanding the regional and species specializations of the human and macaque retina as well as the molecular regulation of aging progression and related diseases.

Single-cell transcriptomes of primate retina
We collected single-cell transcriptomic profiles of 119,520 cells, including 38,558 from six healthy human samples (one 8-day-old infant and, five adults aged between 35 and 87 years) and 80,962 from five macaque samples (one 2-year-old juvenile and four adults aged between 4 and 23 years) (Supplementary Table 1), which covered progressive retinal aging, using droplet-based scRNA-seq (10x Genomics platform). Retinae were dissociated from 1.5-2.0 mm and 5 mm diameter pieces around the center of the foveal pit to represent foveal and peripheral cells, respectively ( Fig. 1a-c, Supplementary Fig. 1a). In total, cell atlases of 45,231 foveal and 74,289 peripheral cells were collected. The sequencing depth was analyzed by the gene number per cell in each sample ( Supplementary Fig. 1b). We performed unbiased clustering of cell profiles from the human retina with both foveal and peripheral regions and found 56 clusters ( Supplementary  Fig. 1c). Cells from the 8-day-old infant generally clustered together and differed from cells from other samples. As human photoreceptors do not reach maturation until 4 years old [19], the cells in the 8-day-old infant were not as mature as those of the adults (Supplementary Fig. 1c). Interestingly, some clusters showed regional specificity, e.g. Clusters 13, 22 and 51 showed foveal identity and Clusters 21, 32 and 43 showed peripheral identity ( Supplementary Fig. 1c). With known cell-type markers, we grouped 56 clusters into 11 major classes, including six neuronal cell types (rods, cones, HCs, ACs, BCs and GCs), three glial cell types (MGs, astrocytes and microglia), endothelial cells and pericytes (Fig. 1b, Supplementary Fig.  1d NSUN4  UNC119  NDUFA4  GUCA1A  NME1  MBD3  GUK1  RPL90  NDUFB9   GOLGB1  AAK1  MAK  MGLL  SCG5  PLCB1  DDX17  RBFOX1  NSUN4  USP15  LIMCH1  Gene loadings   CPE  IFI27L2  RCVRN  NDUFB1  SLIRP  NDUFA1  POLE4  COX7B  UBAP1L  CALR  CASQ1  RS1  CLUL1   Human  Shared  Macaque   15   10   5   0   75   50   25   0   75  50  25   %Expression Expression respectively) showed similar transcriptomic profiles, with high expression of OTX2 and RLBP1 ( Supplementary Fig. 1e). OTX2 is a homeodomain transcription factor expressed in early and mature photoreceptors [5,11], and RLBP1 is a marker for MGs [20]. These two markers have not been found co-expressed in retinal cells previously. To validate the existence of the subtypes, we performed immunofluorescence staining and observed that ∼3% of RLBP1 + cells also expressed OTX2 in the human adult retina, in both foveal and peripheral regions ( Fig. 1d and e), indicating the existence of these cells in the adult human retina. Therefore, Clusters 46 and 51 are previously unidentified retinal cell types.
In the macaque dataset, 68 clusters were classified, and 10 major cell types were annotated (Fig. 1c, Supplementary Fig. 1f and g, Supplementary Table 3). Compared to humans, astrocytes were not detected in the macaque retina (Fig. 1c, Supplementary  Fig. 1f and g), consistent with previous scRNA-seq results from M. fascicularis [3]. Thus, we identified all major retinal cell types and OTX2 + RLBP1 + cells in the primate retina by scRNA-seq analysis.

Cross-species transcriptome comparison between human and macaque retina
To investigate the similarities and differences between the human and macaque retina, we integrated and analyzed the single-cell transcriptomic profiles to perform cross-species comparison using the LIGER algorithm [21], which can identify shared-or specific-'factors,' a cluster of genes that define specific cell types (Fig. 1f, Supplementary Fig. 1h). The major cell types showed high correlations between the two species, suggesting conservativeness between humans and macaques. In addition, GCs and HCs exhibited similar transcriptomes (Fig. 1g), which is possibly a result of their close relationship in development [22]. The LIGER joint clustering assignments indicated that rod clusters showed species-related differences, although they were well-mixed in the tSNE plots ( Fig. 1f and h). Clusters 1, 4 and 15 were humandominant, whereas Clusters 3, 12 and 16 were macaque-dominant (Fig. 1h). We next plotted the species-specific dimorphic genes derived from particular rod clusters (factor 5) in LIGER analysis (Fig. 1i). General rod marker genes (GNGT1 and NRL) were shared, but species-specific genes for rod clusters in humans (CPE, IFI27L2 and NDUFB1) and macaques (GOLGB1, AAK1 and MAK) were also identified (Fig. 1i). Further analysis indicated that these genes also showed cluster specificity, consistent with species preference (Fig. 1j). For example, CPE (carboxypeptidase E), which is required for the maturation of photoreceptor synapse and its signal transmission to the inner retina [23], was highly expressed in human-specific clusters (Clusters 1 and 15) but minorly expressed in some shared clusters (Clusters 10 and 24) (Fig. 1j). To validate the species-specific marker (Fig. 1j), we stained BTG1, a cell cycle inhibitory gene [24], in both human and macaque retina ( Fig. 1k and l). The results showed that BTG1 co-localized with NR2E3 + rods in outer nuclear layer (ONL) retina. However, the ratio of BTG1 − NR2E3 + cells to NR2E3 + cells was much lower in human retina, which is consistent with the scRNA-seq data. Although clusters of the cones and GCs from humans and macaques were closely matched by LIGER assignments (Clusters 23 and 8), we still found several species-specific genes in the cones (factor 18) and GCs (factor 6), respectively ( Supplementary Fig. 1i). Collectively, our data demonstrated transcriptome conservativeness between the two species, but divergent gene expression was also discovered.

Distinct subtypes of primate rod photoreceptors
Rods are the dominant cell type in the mammalian retina and are responsible for vision under low-light intensity [1]. In total, we collected 16,686 rods in human retinal samples. As rods showed the highest interspecies differences by LIGER analysis, we next precisely analyzed rod subtypes. We first removed rod subclusters with more than 98% cells from the 8-day retinal sample to eliminate the influence of development. Based on analysis of differentially expressed genes (DEGs), human rods could be divided into two subgroups based on MYO9A expression ( Fig. 2a and b). CPE, which was found Natl Sci Rev, 2021, Vol. 8, nwaa179 to be dominantly expressed in humans by crossspecies comparison (Fig. 1i), was highly expressed in MYO9A − rods (Fig. 2c). In the human peripheral retina ONL, over 95% of cells were rods. Thus, we performed MYO9A RNAscope in situ hybridization (ISH) of the human peripheral retina and validated that a proportion of human rods expressed MYO9A (Fig. 2d), suggesting that human rods consist of MYO9A + and MYO9A − cells. Consistent with human retina, macaque rods also consist of MYO9A + and MYO9A − cells ( Supplementary Fig. 2b). Moreover, a subcluster of MYO9A − in macaque retina was much lower than that in humans. To investigate the transcriptional regulation of MYO9A, bulk ATAC-seq was performed in human foveal and peripheral retinal samples (27,54 and 86 years old). By searching transcription factor motifs from ATACseq peaks close to the MYO9A transcription start site (TSS), we found two potential OTX2 binding sites in MYO9A gene (Fig. 2e). Consistently, in humans, OTX2 was generally expressed in the MYO9A + subcluster of rods (Fig. 2c, Supplementary Fig. 2a).
Considering OTX2 is a transcription factor controlling retinal photoreceptor cell fate [22], our results indicate that the expression of MYO9A might be regulated by OTX2. Accordingly, we observed both OTX2-positive and -negative cells in the human (Fig. 2f) and macaque retina ( Supplementary  Fig. 2c), similar to the division of MYO9A. From the scRNA-seq data, we found a similar MYO9A + to MYO9A − rod ratio in the foveal and peripheral regions ( Fig. 2g and h). As rods showed species specificity in their transcriptomes, we analyzed the macaque rods (36,904 cells) and observed that one subcluster (4.7% rods) was MYO9A − , which was much lower than that in humans (36.9%) ( Fig. 2i-k) and consistent with ISH data (Supplementary Fig. 2b). Additionally, the distribution of MYO9A − rods showed regional preference, with 98.6% of MYO9A − rods located in the peripheral retina in macaques (Fig. 2l). Taken together, though the rods of humans and macaques could both be divided into two subtypes based on MYO9A expression, the proportion and location of MYO9A − cells differed substantially between the two species. BCs are neurons that connect the outer to the inner layer in the retina, and play direct or indirect roles in transmitting signals from photoreceptors to GCs [1]. We found 17 and 18 BCs subclusters in humans and macaques, respectively, and categorized these subclusters into 12 and 11 known types in humans and macaques, respectively (Supplementary Fig. 2d-f). The 11 types of M. mulatta BCs were consistent with previous findings in M. fascicularis [3]. Similar to M. fascicularis, the Blue Bipolar (BB) and Giant Bipolar (GB) types were grouped together because of their transcriptional similarity in M. mulatta, whereas BB and GB types of human BCs could be separated by the transcriptome with distinct molecular markers ( Supplementary Fig. 2e). However, the BC subtypes were generally conserved between humans and macaques ( Supplementary  Fig. 2g).

Regionally different transcriptomes of primate MGs and cones
The foveal retina is cone-dominated, whereas the peripheral retina is rod-dominated. MGs are particularly important in the fovea because the foveola center is formed only by cone and MG processes. To investigate the region-specific features of MGs and cones, we compared their single-cell transcriptomes in the foveal and peripheral retina. MGs exhibited remarkable regional specificity (Fig. 3a). Generally, most foveal MGs were cells in Clusters 13 and 22 (Fig. 3a) and exhibited high TRH, FGF9, DIO2, CYP26A1, RGR and HTRA1 expression (Fig. 3b,  Supplementary Fig. 3a and b). The majority of MGs in Clusters 21, 32 and 43 were located in peripheral retina and exhibited high NPVF, MT3, MT2A and GPX3 expression (Fig. 3b, Supplementary Fig. 3a  and b). Immunostaining of NPVF and RLBP1 in the human retina validated the scRNA-seq results, indicating that NPVF was a marker gene of peripheral MGs in humans (Fig. 3c). In addition, immunostaining and ISH also revealed the foveal MG specific expression of TRH (thyrotropin releasing hormone) and DIO2 (iodothyronine deiodinase 2), respectively ( Fig. 3d and e). We performed weighted gene co-correlation network analysis (WGCNA) of MG regional DEGs and found several foveal and peripheral MG-specific modules (Supplementary Fig. 3c and d). For the top foveal module, we performed network analysis and found that TRH, RGR, DIO2 and CYP26A1 were also highly enriched (Fig. 3f). CYP26A1 is a retinoic acid-metabolizing enzyme, with foveal-specific expression validated here ( Supplementary Fig. 3e). CYP26A1 has been identified as a DEG of the foveal retina during human development [25], which functions in rod-free zone generation [26,27]. Interestingly, the expression level of RGR (retinal G protein-coupled receptor), another foveal MG-specific gene, is higher in the foveal MGs. RGR encodes a crucial factor for the photic visual cycle of cone visual pigment regeneration, mutations of which are associated with retinitis pigmentosa (RP) [28]. The cone visual cycle of regenerating 11-cis retinal is dependent on both retinal pigment epithelium (RPE) and MGs [28,29]. Thus, the higher RGR expression in the foveal MGs might meet the high demand of cone visual pigment regeneration in the cone-dominated fovea. Additionally, HTRA1 (high-temperature requirement protein A1), a risk factor gene of AMD, exhibited high expression in the fovea MGs, consistent with the fovea being more vulnerable in AMD [30]. Together, our results indicated that MGs in the fovea, as a type of supporting glial cell, expressed specific genes for facilitating cone photoreceptor regeneration, specification and functional maintenance. We next analyzed the MGs of macaques, which also showed regional differences in terms of cell types and gene expression patterns (Supplementary Fig. 3f and g). Comparing the DEGs between the foveal and peripheral MG across both species, several genes showed similar regional specificity, for example TRH, DIO2, CYP26A1, RGR, GPX3 and RPL7, whereas gene expression preference of NPVF, MT1F and JUNB was not conserved (Fig. 3g).
Interestingly, cone photoreceptors also exhibited remarkable region-specificity in both humans and macaques (Fig. 3h). The thyroid hormone (TH) receptor, THRβ (human thyroid hormone receptor beta), is expressed in cones and also serves as a transcription factor [31]. Our ATAC-seq results showed that the potential binding sites of THRβ were found in open chromatin close to the TSS of several foveal-enriched genes (RFC1, POLE4 and  APLP2) (Fig. 3i, Supplementary Fig. 3h). DIO2 converts the inactive TH thyroxine (T4) in blood into the active TH triiodothyronine (T3) [32]. Of note, we also found DIO2 to be highly expressed in foveal MGs (Fig. 3b and e), indicating that foveal MGs may regulate the region-specific transcriptome of cones. The dosage of TH has been reported to play a role in specifying cone subtypes in human retinal organoids [33]. Here, our data demonstrated that MGs might modulate region-specific gene expression by TH signaling in foveal cones in adulthood. CALB1 was one of the markers for peripheral cones, as determined by scRNA-seq data analysis (Fig. 3i). To confirm this finding, we carried out immunostaining of cone OPSIN and CALBINDIN (encoded by gene CALB1) in human retina and found that CALBINDIN was only expressed in cones (both L/M and S cones) located in the peripheral retina (Fig. 3j, Supplementary  Fig. 3i). Taken together, both MGs and cones displayed distinct regional gene expression profiles and MGs may modulate region-specific gene expression to facilitate cone photoreceptor generation, specification and functional maintenance by cell-cell communication.

Molecular classification of retinal HCs in primates
HCs are located in the inner nuclear layer (INL) of the vertebrate retina, where they interconnect laterally with photoreceptors. Here, the human HCs clustered into classical H1 (LHX1 + PCP4 + ) and H2 (ISL1 + CALB1 + ) subtypes ( Fig. 4a and b, Supplementary Fig. 4a). Interestingly, H1 and H2 cell types showed regional preference, i.e. 92.0% of foveal HCs were H1 type and 60% of peripheral HCs were H2 type (Fig. 4c). Immunostaining of CALBINDIN/ONECUT2, and ISH of LHX1 and ISL1 in the human retina confirmed the scRNA-seq results, illustrating that most of the human foveal HCs were H1 (Fig. 4d-f). The macaque HCs were also classified as H1 and H2 types (Fig. 4g, Supplementary Fig. 4b). Compared to the 92.0% in humans, 60.0% of macaque foveal HCs were H1 type (Fig. 4h). Although most DEGs between H1 and H2 showed similar expression across the two primate species, some human H1-enriched genes, e.g. FILIPL1, CMSS1 and PCDH9, were highly expressed in macaque H2 cells, whereas some human H2-enriched genes, e.g. SYNPR, MGARP and AKR1B1, were highly expressed in macaque H1 cells (Fig. 4i, Supplementary Fig. 4c), indicating that distribution and gene expression were slightly different across the H1/H2 cell types of both species.

Aging patterns of human retina
Aging affects retinal function and structure physiologically and pathologically. To clarify the cellular and molecular changes during aging, we systematically studied human retinal samples (between 35 and 87 years old). We categorized the five samples into three groups (35Y, adult group; 52Y and 63Y, mid-age group; 86Y and 87Y, aging group). We first focused on whether the genes continuously increased or decreased from the adult to aging groups (Fig. 5a-e). In total, the expression of 87 genes was up-regulated. Gene ontology (GO) analysis indicated these genes may play roles in response to hypoxia, regulation of cell death, microglial cell activation ( Fig. 5a and b, and Supplementary  Table 4). Interestingly, some up-regulated genes only increased in certain cell types. For example, LIN7A and CALD1 were only elevated in cones with aging (Fig. 5c). Additionally, 121 genes related to visual perception, phototransduction, adenosine triphosphate (ATP) biosynthetic process and retinoid metabolic processes, were down-regulated with aging ( Fig. 5d and e, Supplementary Table 4), indicating that visual function is affected as retina ages. We also analyzed the age-related GO change in each cell subtype ( Supplementary Fig. 5a-h). Interestingly, in the GO analysis for each cell type, we found that most retinal neurons showed similar GOs to total retinal cells (Fig. 5b), such as response to hypoxia, response to oxidative stress and regulation of cell death. In contrast, GOs in MG are antiapoptosis, negative regulation of programmed cell death and response to hypoxia. These data suggest that retinal cells are heterogeneous during aging. The GO analysis suggested that MG might be in antiapoptosis status to prevent retinal neuron death during aging. Consistently, among the neurons, a dramatic decrease in rod photoreceptors with retinal aging was observed by scRNA-seq data and immunostaining analysis ( Fig. 5f and g), but there was no dramatic cell number change of BCs (Supplementary Fig. 5i and j). Rod outer segments, which are specialized compartments of photoreceptor for phototransduction, were greatly decreased with RHO protein (Fig. 5g). In addition, RHO protein was ectopically expressed in cell bodies, indicating defective light-sensing function of aged rod photoreceptors. In retinal glial cells, the proportion of microglia was remarkably increased in the aged samples ( Fig. 5f and h), consistent with the upregulated genes, which were enriched in response to hypoxia and microglial cell activation (Fig. 5b). Interestingly, the major changes in age-related cell-type proportion in humans were not obvious in the macaque retina (Fig. 5f).  As two major classes of rods were identified in humans (MYO9A + and MYO9A − ) (Fig. 2) and rods were the most vulnerable cell type during aging, we next asked whether these two types of rods were differentially affected by aging. Quantification analysis of scRNA-seq data indicated that the population of MYO9A − rods dramatically decreased during aging (Fig. 5i), which was confirmed by the ISH of MYO9A in the human peripheral retina ONL (Fig. 5j). To further validate the reduction of MYO9A − rods, we quantified the ratio of MYO9A + and MYO9A − cells to the total retinal cells in different ages of human retinae (Supplementary Fig. 5k). The data revealed that both MYO9A + and MYO9A − cells were reduced along aging. Moreover, the proportion of MYO9A − cells was much more dramatically reduced in aged human retina than was MYO9A + cells. Thus, it appeared that human MYO9A + rods were more resistant to aging. Because the macaque rods were predominately MYO9A + (95.4%), we did not observe obvious changes in the proportions of MYO9A + and MYO9A − rods in macaque retina during aging (Fig. 5i). As OTX2 potentially regulates MYO9A expression in human rods (Fig. 2e), we next explored OTX2 expression in the human peripheral retina ONL. Unsurprisingly, OTX2 + cells in the peripheral ONL increased from 62.1% at 54 years old to 96.7% at 87 years old (Fig. 5k and l). Together, we identified a cell-type-specific degeneration of

Aging-related changes in cell-cell communication
Aging is a systematic continuous process, regulated by both intrinsic and external signals. To understand how intercellular communication plays a role in the aging of rods, we built a network of ligand-receptor interactions among different cell types and focused on the changes in ligand-receptor expression with aging (Fig. 5m, Supplementary Fig. 5l and m). We studied the ligand-receptor changes between rods and MGs or BCs (Fig. 5m, Supplementary  Fig. 5m), as MGs and BCs closely interact with rods. Expression of VEGF-A, which are linked to AMD and retinal angiogenesis, were up-regulated in BCs and MGs. Their receptors also increased in rods with aging (Fig. 5m). Interestingly, VEGF-A expression was dramatically increased in MGs, HCs and astrocytes. Accordingly, abnormal angiogenesis was also observed in the aged retina ( Fig. 5n and o,  Supplementary Fig. 5n). The length of blood vessels was significantly increased in aged retina. Together, these findings indicate that aging could be a synergistic process with comprehensive changes in different types of neuronal or glial cells in the human retina.

Molecular changes in foveal and peripheral aged human retina
To investigate the aging trajectories of foveal and peripheral primate retina, we established aging scores using an aging-related gene set which was downloaded from Human Ageing Genomic Resources (https://genomics.senescence.info/). The methods used for aging score calculation were adapted from the algorithm used by Nowakowski to predict the maturation of cortical neurons [34] (see Methods). We found that the elevation trends of the aging score generally matched the real sample age (Fig. 6a). Of note, the foveal cells showed higher aging scores than peripheral cells, indicating that the fovea may be more vulnerable to aging. In addition, the increase in aging scores started after 52 years old (Fig. 6a), suggesting that the retinal aging process may gradually accelerate at middle age. To verify the reliability of the aging score, we first examined the expression trends of classical aging genes, and found that the expression levels of STAT3 and PIK3R1 were positively correlated with aging. However, the expression levels of PCNA and GSTP1 were negatively correlated with aging, indicating that the aging score models by single-cell transcriptomic profiles were reliable (Fig. 6b, Supplementary  Fig. 6a). We then focused on the expression of human aging up-regulated genes in the retina and found that NPVF, GPX3, B2M, PCP4 and CRYAA were positively correlated with aging, whereas down-regulated genes RHO, DRD4, GNGT1 and MFG8 were negatively correlated with aging (Fig. 6c, Supplementary Fig. 6b and c, Supplementary Table 4).
As foveal and peripheral retina exhibited different degrees of aging (Fig. 6a), and MGs, cones, and HCs showed regional differences (see Figs 3 and 4), we investigated whether these regional differences changed during aging. NPVF, a highly expressed gene in peripheral MGs, was significantly up-regulated during aging based on the scRNA-seq c d e f Figure 7. Region-, cell type-, and species-specific expression patterns of genes associated with human retinal diseases. (a and b) Aggregated expression of human eye disease-associated genes in human (a) and macaque (b) cell types, with some diseases showing obvious regional enrichment characteristics. (c and d) Aggregated expression of age-related macular degeneration-associated genes at different human (c) and macaque (d) ages.
(e and f) Expression patterns of human eye disease-associated genes in human (e) and macaque (f) retinal cell types. P values were calculated by bootstrap hypothesis test, * P < 0.05, * * P < 0.01. dataset (Fig. 6d). Consistently, immunostaining of NPVF suggested its increased expression in peripheral MGs with aging. Interestingly, NPVF was ectopically expressed in the foveal MGs of 86-yearold retina, but not of 54-year-old retina (Fig. 6e). Additionally, we found that average expression of foveal MG-enriched genes gradually decreased during aging, indicating that the regional specificity of foveal MGs declined ( Supplementary Fig. 6d). In contrast, the expression of peripheral MG-enriched genes did not show an obvious change trend during aging ( Supplementary Fig. 6e).
Other than MGs, the proportion of H2-type HCs in humans decreased from 52.3% (mid-age) to 14.2% (aged) (Fig. 6f). However, this trend was not observed in macaques ( Supplementary Fig. 6f), indicating that H2 cells were more sensitive to aging processes in humans than in macaques. We also used immunostaining to verify the changes in HC cell type with aging. The number of H2-type cells decreased in both regions, whereas H1-type cells showed no change ( Fig. 6g and h). By comparing the DEGs of H1 and H2 HCs, we found that HSP90AA1 and HSP90AB1, two heat shock proteins, were highly expressed in H1 but not H2 (Fig. 6i). These proteins can mediate lysosomes to refold unfolded proteins, thus protecting cells from the effects of protein toxicity [35]. This may be a reason why the H1 cells showed higher anti-aging ability than H2 cells.
For cones, the regional expression and localization of CALBINDIN showed no differences with aging. However, the subcellular localization of cone OPSINs changed from the outer segments to cell bodies in aged cones (both L/M and S cones) (Fig. 6j-m). It is possible that dysfunction of protein transportation in cilia leads to this mistrafficking and accumulation of cone opsins [36].

Human retinal disease related to gene expression in retinal cells
We used our scRNA-seq dataset to investigate the cell-type and region specificity of 178 genes associated with 55 types of human retinal disease, including night blindness, macula dystrophy, rod-cone dystrophy, dominant RP, recessive RP, AMD and recessive achromatopsia (RETNET, https://sph.uth.edu/retnet/). We first calculated the expression of each gene in the foveal or peripheral cell clusters and then aggregated the expression scores by disease types (Fig. 7a and b,  Supplementary Fig. 7a and b). Generally, rods and cones, especially foveal cones, were highly related to night blindness and macula dystrophy ( Fig. 7a and  b). Genes related to rod-cone dystrophy, and dominant and recessive RP were highly expressed in rods and cones, particularly foveal photoreceptors in both humans and macaques. Genes related to AMD were highly enriched in foveal MG and cones, suggesting a correlation of regional cell subtype with this disease.
We further illustrated cell-type and region specificity of each gene categorized by disease. As expected, age-related degeneration AMD-risk genes were highly expressed in aged human retina but not in that of aged macaques ( Fig. 7c and d), suggesting fewer aging characteristics in the macaque retina. Cell-type and region specificity of disease-related genes were further compared between humans and macaques ( Fig. 7e and f, Supplementary Fig. 7a and b). For example, PDE6A and PDE6B, which are reported to cause recessive RP, were conservatively expressed in the rods in both primates ( Fig. 7e and f). HTRA1 which is a risk factor for AMD, were highly expressed in HCs in both humans and macaques but with differential regional preference ( Fig. 7e and f). Another recessive RP-associated gene, MVK (mevalonate kinase), was enriched specifically in human peripheral HCs ( Fig. 7e and f). Mutations of EFEMP1 can lead to macular dystrophy [37], and the expression of EFEMP1 was high in human and macaque foveal MGs (Fig. 7e and f). In addition, ATF6, a recessive achromatopsia-associated gene [38], was selectively expressed in human GCs located in the periphery, but showed no obvious expression in macaques ( Fig. 7e and f). Together, these datasets showing region-and cell-type-specific expression patterns of retinal-disease-related genes should be useful for mechanistic and therapeutic studies of retinal diseases in the future.

DISCUSSION
By profiling 119,520 single cells, we provided detailed single-cell transcriptomes of adult and aged primate retinae with regional information. The retinal transcriptome at single-cell or single-nuclei level has been reported [3,[7][8][9][10][11][12]. In comparison with other single cell retina studies, we have the largest number of primate retinal samples in one study (six human samples and five macaque samples) covering the longest life span (8 days after birth to 87 years old for humans; 2-23 years old for macaques) so far. The non-human primate we used here is Macaca mulatta, which is evolutionarily closer to humans than Macaca fascicularis [3], adding more primate retinal transcriptome data to the retinal resource. Moreover, we did cross-species analysis between humans and macaques, and found the species-specific signatures of primate retinae, which provided valuable data sources to understand the regional and species specializations of the human and macaque retina. In particular, the study provided conceptual advances in the molecular characteristics of aging progress of human retina. Our data revealed that human retinal aging occurred in a region-and cell-type-specific manner, suggesting strong cell heterogeneity in retinal aging. Moreover, a big dataset was generated and showed the cell-type-, region-, and speciesspecific gene expression associated with various human retinal diseases, which provides a foundation to understanding of the molecular and cellular mechanisms underlying human retinal diseases. Overall, this study provides a valuable resource at the single-cell resolution to help understand human and macaque retina and related aging processes, study the regulatory mechanism of the aging process of human retina, and identify the molecular markers for aging and degenerative changes of human retina.

Different rod subtypes in human and macaque retina
We examined retinal cell-type evolutionary conservation between humans and macaques. The major cell types were conserved in the two species, but the molecular features in retinal cells showed species differences (Fig. 1f-k). Based on the expression level of MYO9A, the rods in humans and macaques can be divided into two groups. Our data indicated that the MYO9A − rod subtype was reduced during aging ( Fig. 5i and j), which may result from either depletion of MYO9A − rods or increased MYO9A + cells. Considering that the proportion of MYO9A − cells was much more dramatically reduced in aged human retina than was the MYO9A + cells (Fig. 5i,  Supplementary Fig. 5k), it is highly possible that the reduction of rods was largely a result of depletion of MYO9A − rods. Future studies are needed to investigate the functional role of MYO9A in rod photoreceptors. It is worth noting that we identified previously unidentified cell types. Cells in Cluster 46 of rod identity and Cluster 51 of MG identity co-expressed OTX2 and RLBP1 (Fig. 1d and e,  Supplementary Fig. 1e). OTX2 is a typical marker for BCs and rods, and RLBP1 is a marker for MGs. MGs are considered to be retinal precursor cells in non-mammalian vertebrates and are able to transdifferentiate into rods and BCs with additional transcriptional and epigenetic regulation in mice [39,40]. The co-expression of OTX2 and RLBP1 suggests that MGs in Cluster 51 may be transcriptionally closer to rods and BCs. However, their function and potential to transdifferentiate into neuronal cells, such as rods or BCs, need further exploration.

Distinct transcriptomes of primate foveal and peripheral retina
With enriched cones, the fovea is the most important region for high-acuity vision in humans, whereas the peripheral retina helps direct eye movements to focus salient images on the fovea [41]. Based on their different functions, they showed distinct molecular signatures in MGs and cones (Fig. 3a, b, h and i). Several genes were found to be specifically enriched in foveal MGs, for example RGR, HTR1, DIO2 and CYP26A1. The finding of specific high RGR expression in foveal MGs is interesting. RGR is a non-visual opsin in intracellular membranes of RPE and MGs. Previous study suggests that RGR opsin and retinol dehydrogenase-10 (Rdh10) convert alltrans-retinol to 11-cis-retinol for the regeneration of cone visual pigment in daylight, as RPE cannot meet the higher needs of the visual opsins in rods and cones in daylight [28]. An intriguing speculation is that the higher expression of RGR in fovea MGs is initiated to meet the higher demands of visual pigment regeneration in foveal cones.
In addition to RGR and retinol, we also found that MGs and cones possibly interact with each other via the TH signaling pathway. As the active form of TH, T3 is important for cone development in retina and retinal organoids [33,42]. However, it remains unclear whether TH signals exhibit regionally specific distribution in primate retina, and which cell types may provide TH. The major TH in blood is inactive T4, which needs to be converted to T3 by DIO2. Interestingly, we found DIO2 to be selectively expressed in foveal MGs. Accordingly, we found that cone photoreceptors expressed THRβ, a nuclear receptor of T3. Based on our ATAC-seq data, we observed the binding sites of THRβ in the promoter regions of RFC1, APLP2 and POLE4, which are genes specifically expressed in foveal cones. Hence, it is possible that the expression of DIO2 in foveal MGs may produce high local T3, which modulates the expression of a specific set of genes in foveal cones. Of note, higher expression of NPVF, which is negatively regulated by the TH signaling pathway [43], was detected in the peripheral cones and MGs. Together, our data imply that foveal MGs may provide a local TH signal that regulates the cellular specificity of foveal and peripheral cones.

Cellular and molecular changes accompanying human retinal aging
Our study explored cellular and molecular changes during primate retina aging. GO term analysis revealed that microglial cell activation, protein refolding, oxygen species metabolism and cell death were highly enriched in aging retina and could possibly be considered as hallmarks for retinal aging. In healthy retinae, microglial cells are innate immune cells, which can constantly adapt to their microenvironment. In mouse models, increased density and activation of microglia exist in aged retinae and brains, as well as in age-related human diseases, including AMD, Alzheimer's disease (AD) and Parkinson's disease (PD) [44,45]. Similarly, the human retina showed an increased number of microglia during aging (Fig. 5h). Moreover, microglial activation and microglia-mediated inflammatory responses can generate a chronic mild inflammatory environment, including an increased production of inflammatory cytokines, reactive oxygen species (ROS) and reactive nitrogen species (RNS) [46], which can further accelerate retinal aging and pathogenesis of age-related diseases. Thus, we speculated that targeting microglial activation states could potentially be used as a therapeutic avenue in slowing aging and retinal diseases.
Mistrafficking and accumulation of both rod and cone opsins (Fig. 5g and Fig. 6j and l) were found in the aging retina. Damaged protein accumulation can lead to toxicity against photoreceptors as these failed homeostatic mechanisms can contribute to pathology [36]. Therefore, protein misfolding may be a causative reason for apoptosis of photoreceptors during aging. Other than protein misfolding, retinal oxygen supply also decreases with age. Thus, aging retinae suffer from low-grade chronic ischemia [47], which can lead to retinal oxidative stress in turn. MGs in the peripheral retina exhibited high expression of metallothioneins (MT3, MT1G and MT2A) and glutathione peroxidase 3 (GPX3), which participate in an array of protective stress responses [48]. Metallothioneins are a family of cysteine-rich metal-binding proteins involved in protection against oxidative stress and buffering against toxic heavy metals [49]. As a major scavenger of ROS in plasma, GPX3 acts as a redox signal modulator [50]. The high levels of metallothioneins and GPX3 in the peripheral MGs may contribute to fovea to peripheral aging gradient observed in human retina.
VEGF is a well-known risk factor for abnormal retinal angiogenesis in AMD and diabetic retinopathy. We found that VEGF-A expression increased specifically in HCs, MGs and astrocytes following aging. VEGF is a potent endothelial cell mitogen that stimulates proliferation, migration and tube formation, leading to angiogenic growth of blood vessels [51]. With the knowledge that cellular origin of VEGF-A increases during aging, it is possible to design therapeutic approaches to reduce its expression in MGs, HCs and astrocytes for abnormal retinal angiogenesis treatment. Other than VEGF, we observed several previously unidentified genes associated with aging, for example NPVF, PCP4, CRYAA and B2M. Altogether, our study provides a valuable data source for future studies on human retinal aging and related diseases.

Cell-type and region specificity of human retinal diseases
Our study provides a foundation to understand the molecular and cellular mechanisms underlying human retinal diseases. By assessing the expression of 178 genes implicated in human retinal diseases, we found that these genes were preferentially expressed in specific retinal cell classes and regions. Photoreceptors in foveal regions are the major cell types for many retinal degenerative diseases. These findings emphasized the value of human aging transcriptomes for age-related retinal diseases. In particular, AMD-risk genes were highly expressed in the aged human retina. Together, the regionand cell-type-specific expression patterns of retinal disease-related genes may provide guidelines for future cellular and molecular precision therapies.

Human retinal samples
The isolation procedure of human retinal samples and research protocols were approved by the Research Ethics committee of the Peking University Third Hospital and Peking Union Medical College Hospital, and were conducted in accordance with approved institutional guidelines. All the protocols are in compliance with the 'Interim Measures for the Administration of Human Genetic Resources,' administered by the Ministry of Science and Technology of China.

Macaque retinal samples
Rhesus macaques (Macaca mulatta) were provided by Xieerxin Biology Resource (Beijing, China). All procedures were conducted in accordance with the Principles for the Ethical Treatment of Non-Human Primates and were approved by the Institutional Animal Care and Use Committee of the Institute of Biophysics and University of Science and Technology of China. All animals were maintained at 25 • C on a 12 hour light:12 hour dark schedule at Xieerxin Biology Resource, a Laboratory Animal Care-accredited facility in Beijing, in compliance with all local and federal laws governing animal research. Animals were given a commercial diet twice a day with tap water provided ad libitum and were fed vegetables and fruits once daily under careful veterinary oversight. Before the experiment, none of the animals had a clinical or experimental history that would affect physiological aging or increase disease susceptibility.

Statistical analyses
All data were represented as means ± s.e.m. Comparisons between two groups were made using t-tests. Quantification graphs were analyzed using GraphPad Prism (GraphPad Software).

DATA AND CODE AVAILABILITY
The single-cell RNA-seq and bulk ATAC-seq data used in this study were all deposited in the Genome Sequence Archive (GSA) for Human, National Genomics Data Center (https://bigd. big.ac.cn/gsa-human) under BioProject accession numbers PRJCA002731. The GSA number for macaque is CRA002680 and the GSA number for Human is HRA000182. All codes for data analysis are included in Supplementary file 1.