Preparation of human neural retina samples and generation of single‐cell transcriptome atlas
We obtained post‐mortem human adult eyes approved for research purposes following corneal transplantation. As the transcriptome profile of human retinal pigment epithelial cells has already been reported (Liao
et al,
2010; Strunnikova
et al,
2010), we focused solely on the neural retina layers. In this study, we extracted the neural retina from 12 donor eyes (
Appendix Table S1). We observed consistent cell viability across retinal tissues retrieved within 15 h post‐mortem (
Appendix Fig S1A) and found that donor age does not impact negatively on cell viability in the extracted neural retina (
Appendix Fig S1B). To minimize potential risk of mRNA degradation due to reduced cell viability, we selected three donor samples retrieved within 15 h post‐mortem and analyzed them with single‐cell RNA sequencing (scRNA‐seq) using the 10X Genomics Chromium platform.
Sequence data from five scRNA‐seq libraries derived from the three neural retinal samples were pooled for processing and analysis. From 23,000 cells, we obtained an average of 40,232 reads per cell and 1,665 UMIs (unique transcripts) per cell. Following quality control and filtering using the Seurat package (Butler
et al,
2018), our final dataset contained 20,009 cells, which were taken forward for further analysis.
The scRNA‐seq data were initially analyzed using an unsupervised graph clustering approach implemented in Seurat (version 2.2.1) to classify individual cells into cell populations according to similarities in their transcriptome profiles. Overall, the cells were classified into 18 transcriptionally distinct clusters (
Appendix Fig S2). We first assessed the variation between donor samples (
Appendix Table S2). Interestingly, although many of the identified clusters are well represented in all three donor retinal samples, we also observed several donor‐specific clusters corresponding to rod photoreceptors (
Appendix Fig S3A). In contrast, we observed minimal variation between two different libraries prepared from the same donor sample, supporting the quality of the scRNA‐seq datasets in this study (
Appendix Fig S3B). The average expression for all detected genes in each cluster is listed in
Dataset EV1.
Identification of major cell types in the human retina using scRNA‐seq
Based on known markers (Blackshaw
et al,
2001; Imanishi
et al,
2002; Corbo
et al,
2007; Soto
et al,
2008; Klimova
et al,
2015; Macosko
et al,
2015; Shekhar
et al,
2016; Vecino
et al,
2016), we were able to assign cell identities to 16 of the 18 clusters (Fig
1A–D), corresponding to rod photoreceptors (
PDE6A, CNGA1, RHO), cone photoreceptors (
ARR3, GNGT2, GUCA1C), Müller glia (
RLBP1/CRALBP), retinal astrocytes (
GFAP), microglia (
HLA‐DPA1, HLA‐DPB1, HLA‐DRA), bipolar cells (
VSX2, OTX2), retinal ganglion cells (
NEFL, GAP43, SNCG), amacrine cells (
GAD1, CALB1, CHAT), and horizontal cells (
ONECUT1, ONECUT2). The expression of selected marker genes is displayed in
t‐SNE plots (Fig
1D). Two clusters (C5 and C14) express markers from multiple retinal cell types (
Appendix Fig S4); thus, we were unable to assign cell identities to these two clusters and they were excluded from further analysis. Interestingly, our data demonstrated multiple transcriptionally distinct clusters within the rod photoreceptors (six clusters) and bipolar cells (three clusters). In contrast, only one cluster was detected for cone photoreceptors, Müller glia, retinal ganglion cells, horizontal cells, amacrine cells, retinal astrocytes, and microglia, respectively. Correlation analysis confirmed the similarity between clusters within the same cell type (Fig
1E). As expected, we observed high correlations between the expression levels of transcripts within photoreceptor cell types (rod and cones), as well as glial cells (retinal astrocytes and Müller glia) and other retinal neurons (bipolar cells, retinal ganglion cells, amacrine cells, and horizontal cells). The composition of cell populations across our three donors shows that the majority of the cells in human neural retina were rod photoreceptors (~74%) followed by bipolar cells (~10%). These results are similar to those reported in mice, where rod photoreceptors and bipolar cells form the majority of cells in the retina (Jeon
et al,
1998; Macosko
et al,
2015).
To identify genes whose expression was specific to a given cell type, we performed differential gene expression analysis to identify marker genes for each cluster (Fig
1F). We subsequently extracted membrane‐related proteins from gene ontology annotations to identify potential surface markers, which can be used to develop immuno‐based methods to isolate primary culture of individual retinal cell types.
Appendix Table S3 lists the identified markers for individual retinal cell types. We also assessed the gene expression of a panel of commonly known markers in amacrine cells and bipolar cells (Figs
EV1 and
EV2), as well as a panel of markers for subtype identification recently identified in mouse scRNA‐seq studies (Macosko
et al,
2015; Shekhar
et al,
2016). In particular, the bipolar clusters can be classified as OFF‐bipolar cells (
GRIK1+: C6) and ON‐bipolar cells (
ISL1+: C8, C11). Further analysis showed that C8 represents rod bipolar cells based on the marker
PRKCA, while C11 expresses the marker
TTR corresponding to a diffuse bipolar subtype DB4 (Fig
EV1). In summary, we profiled the transcriptomes of all major cell types in the human retina in the presented dataset. Due to their abundance, for the subsequent analyses we focused on the photoreceptors and glial cells.
Profiling healthy and degenerating human rod photoreceptor subpopulations
We profiled 14,759 rod photoreceptors and showed that they can be classified into six populations with distinct gene expressions (C0, C1, C2, C3, C4, and C7). We assessed these six clusters with a panel of seven known rod or pan‐photoreceptor markers (Fig
2A). Our results suggest differential expression patterns among the seven markers. All seven rod markers are highly abundant, consistent with previous scRNA‐seq studies of mouse and human retina (Macosko
et al,
2015; Phillips
et al,
2018). The seven markers showed differential expression patterns in the six identified rod photoreceptor clusters. In particular,
RHO, GNGT1, and
SAG have the highest levels of rod marker detected, followed by
NRL,
ROM1,
GNAT1, and
CNGA1. We also noted that
ROM1 is expressed in both rod and cone photoreceptors, which is consistent with previous studies (Boon
et al,
2008). Importantly, many rod photoreceptor clusters consist of a majority of cells from a single donor (> 90% for C0, C2, and C4 and > 80% for C1 and C7; Fig
2B). It is possible that this observation is due to the systematic biases such as differences in tissue retrieval time, age of donors, or other sample preparation variation. The exception is cluster C3, which is well represented by all three donors.
Next, we set out to further define and classify heterogeneity in rod photoreceptors. We observed that
MALAT1, a long non‐coding RNA that plays a role in retinal homeostasis and disease (Wan
et al,
2017), was robustly expressed in ~66% of the identified rod photoreceptors (9,750 cells), while the rest had little to no expression (5,009 cells; Fig
2C). As such, we utilized
MALAT1 expression as a discriminator and investigated differences between rod photoreceptors with high expression (
MALAT1‐
hi; > 4.68 normalized transcripts per cell) or low expression (
MALAT1‐lo; < 4.68 normalized transcripts per cell).
MALAT1‐
hi and
MALAT1‐lo rod photoreceptors were consistently found in each donor and library samples, with
MALAT1‐
hi accounting for ~66, 90, and 36% of the rods in donors #1, #2, and #3, respectively (Fig
2D). To further validate this finding, we performed RNA
in situ hybridization in another three donor retinal samples. We consistently observed the presence of
MALAT1‐hi and
MALAT1‐lo subpopulations of rod photoreceptors in all retinal samples (Figs
2E and
EV3A). Together, these results showed the presence of heterogeneity within rod photoreceptors that can be distinguished by
MALAT1 expression.
To rule out the possibility that the presence of
MALAT1 rod subpopulations was due to donor sample variations, we applied canonical correlation analysis (CCA) to correct for technical and batch artifacts. We found that CCA effectively corrected the donor‐specific effect on rod photoreceptor clusters (Fig
3A and B). The average expression for all detected genes in each cluster is listed in
Dataset EV2. Following CCA correction, we identified three rod photoreceptor clusters (CCA0, CCA1, and CCA10), which expressed a panel of seven rod photoreceptor markers and were well represented in all donor samples (Fig
3C). Notably, the majority of cells in CCA0 were low in
MALAT1 expression, while CCA1 and CCA10 represented
MALAT1‐hi rod subpopulations (Fig
3D and E). This is consistent with our RNA
in situ hybridization analysis, where we consistently observed
MALAT1‐hi and
MALAT1‐lo subpopulations of rod photoreceptors in all retinal samples (Fig
EV3A). Collectively, these results provide evidence that
MALAT1 heterogeneity in rod photoreceptors is not due to inter‐individual variability.
We also considered the possibility that
MALAT1‐lo rod subpopulations may represent an artifact of “low‐quality cells” in scRNA‐seq data, due to a low number of sequencing reads or broken cell membrane. In this regard, upregulated levels of mitochondrial‐encoded genes and ribosomal proteins can be used to identify such low‐quality cells in scRNA‐seq data (Ilicic
et al,
2016). For our scRNA‐seq dataset, we did not observe upregulation in gene expression for a panel of ribosomal proteins (
RPL41, RPLP1, RPL21, RPS27, RPL13A, RPL36, RPL39, and
RPS28; Fig
3F). However, the rod cluster CCA10, representing 1.4% of rod photoreceptor cells, showed markedly increased levels of mitochondrial‐encoded genes (
MT‐CO2, MT‐ND5, MT‐ND3, MT‐CYB, MT‐ND1, MT‐ND2, MT‐CO3, MT‐ATP6, MT‐CO1, and
MT‐ND4; Fig
3G), suggesting that CCA10 represented a low‐quality
MALAT1‐hi rod cluster and was excluded from further analysis.
As we utilized post‐mortem retinal samples in this study, we reasoned that
MALAT1‐lo subpopulations may potentially reflect the early stages of post‐mortem degeneration in rod photoreceptors. To determine this, we extracted retinal samples from the same donor at different time points of progressive post‐mortem degeneration, with longer time points predicted to have more stressed/dying photoreceptors. Our results showed that there were a high proportion of
MALAT1‐hi rod photoreceptors at 7 h post‐mortem (Fig
4, ~95%). However, we observed a marked decrease in
MALAT1 expression in rod photoreceptors at 13 h post‐mortem. Similar results were observed for the three retinal samples processed for scRNA‐seq (Fig
EV3B). Together, these results demonstrated that
MALAT1 is a novel marker for healthy photoreceptors with loss of expression potentially preceding putative cell degeneration. In summary, we showed that scRNA‐seq can be used to profile healthy (CCA1) and degenerating rod photoreceptors (CCA0), which can be distinguished by high or low
MALAT1 expression levels, respectively.
Transcriptome profile of cone subtypes in the human retina
We detected 564 cone photoreceptor cells in our scRNA‐seq data, which are distinguishable from the other cell types by the expression of the cone marker genes
ARR3, CNBG3, GNAT2, GNGT2, GRK7, GUCA1C, PDE6C, PDE6H, OPN1LW, RXRG, and
THRB (Fig
5A). All 11 marker genes analyzed show specific expression patterns in the cone cluster (C10). We set out to further assess the composition of the cone cluster. In humans, there are three identified subtypes of cone photoreceptor, which can be distinguished by expression of a sole opsin gene:
OPN1SW‐positive S‐cones,
OPN1MW‐positive M‐cones, and
OPN1LW‐positive L‐cones respond preferentially to shorter, medium, and longer wavelengths responsible for color vision (Roorda & Williams,
1999). Notably,
OPN1LW and
OPN1MW exhibit ~98% sequence homology and are unable to be distinguished by 3′ sequencing utilized in this study. By quantifying the number of cells that express the opsin genes, our results showed that the majority of the cone cluster are L/M‐cones (73.22%) and S‐cones in much lower proportion (3.19%; Fig
5B), at levels consistent with those estimated by a previous adaptive optics and photobleaching study (Roorda & Williams,
1999). As expected, the identified cone photoreceptors only express either
OPN1SW or
OPN1LW/MW (Fig
3C).
To further study the molecular differences and identify molecular markers for the cone subtypes, we performed differential gene expression analysis to determine genes that can distinguish the cone subtypes. Our results identified a panel of genes that differentially marked S‐cones (e.g.,
CCDC136,
RAMP1,
LY75, CADM3, TFPI2, CRHBP, RAB17, UPB1, RRAD, and
SLC12A1) and L/M‐cones (e.g.,
THRB, KIF2A, LBH, PGP, CHRNA3, AHI1, and
LIMA1; Fig
3D). We compared this list of cone subtype genes to those identified in scRNA‐seq studies of the macaque and mouse retina, and showed that a number of the cone subtype genes in humans are conserved in macaque and mouse (Macosko
et al,
2015; Peng
et al,
2019), including S‐cone marker
CCDC136 and L/M‐cone marker
THRB. Interestingly,
CCDC136 is located next to the
OPN1SW locus in humans and could possibly be co‐regulated at the transcriptional level. The thyroid hormone receptor
THRB is required for the development of M‐cones in mice (Ng
et al,
2001) and L/M‐cones in humans as determined by a pluripotent stem cell model (Eldred
et al,
2018). Notably, there are two known receptor isoforms for THRB (TRβ1 and TRβ2) and further research to determine the precise roles of THRB isoforms in subtype specification of human cones would be of great interest. Moreover, the transcription factor
TBX2 has been implicated in subtype specification of
Sws1‐cones in zebrafish and chicken (Alvarez‐Delfin
et al,
2009; Enright
et al,
2015). In support of these studies, our data showed that
TBX2 marks the S‐cones in humans, which is also conserved in macaque (Peng
et al,
2019). Together, these results detailed the molecular profiles and identified marker genes that can distinguish the cone subtypes in humans.
Assessment of glial cells in human retina
Next, we focused on two related glial cell types in the human retina, the Müller glia, and the retinal astrocytes. Our scRNA‐seq data have profiled a total of 2,723 Müller glial cells, which classified into a single cluster (C9), and 49 retinal astrocytes, which form a single cluster (C16). Figure
5E shows the expression of a panel of 9 commonly used markers for Müller glia and retinal astrocytes. Our results demonstrated that many of these markers are present in both Müller glia and retinal astrocytes at differential expression levels.
VIM, GLUL, and
S100A1 marked both Müller glia and retinal astrocytes at high expression levels.
GFAP represents a reliable marker for retinal astrocytes, and its expression is consistent with a previous report (Vecino
et al,
2016). Notably, Müller glia are low in
GFAP expression, indicating they are not in an activated state commonly caused by stresses and reactive gliosis (Fernández‐Sánchez
et al,
2015). The
S100B is also expressed in retinal astrocytes at varying levels but absent in Müller glia. Conversely, Müller glia can be distinguished from retinal astrocytes by high expression levels of
RLBP1, and
RGR to a lesser extent. Together, these results provide insights into the differential expression patterns of known glial markers in Müller glia and retinal astrocytes in humans.
As glial cell proliferation has been linked to a range of pathological conditions including retinal gliosis and retinal injury (Subirada
et al,
2018), this provides a means of assessing the health of the profiled retinas. We assigned a cell cycle phase score to each cell using gene expression signatures for the G1, S, G2, and G2/M phases (Kowalczyk
et al,
2015; Fig
EV4). We determined that most of the Müller glial cells expressed genes indicative of cells in G1 phase (75%), suggesting they are not proliferative. These results demonstrated the absence of hallmarks of gliosis/retinal injury and support the quality of the donor retinas profiled.
Using the human neural retina transcriptome atlas for benchmarking
To demonstrate the use of our dataset as a benchmarking reference, we compared the scRNA‐seq profiles of distinct cell types generated using alternative methods, including fetal human cone photoreceptors, human‐induced pluripotent stem cell‐derived cone photoreceptors (hiPSC‐cone; Welby
et al,
2017), and a sample of adult human retina with 139 cells (Phillips
et al,
2018). Correlation analysis demonstrated that the adult human retina sample showed highest similarity to rod photoreceptor (0.63; Fig
EV5), which is expected as rod photoreceptors represent the majority of the cells in the retina. Interestingly, our results also showed that the transcriptome of hiPSC‐cone after 15 weeks of differentiation exhibited the highest similarity to cone photoreceptors, and low similarities to all other retinal cell types (Figs
6A and
EV5). In particular, hiPSC‐cone showed high similarities to fetal cone photoreceptors and adult cone photoreceptors (0.71 and 0.61, respectively), and a much lower similarity to adult rod photoreceptors (0.33). In support of this, principal component analysis demonstrated that the hiPSC‐cone are closer to fetal cone photoreceptors, rather than the adult counterpart (Fig
6B). These results confirmed direct differentiation of hiPSCs to cone photoreceptors with good quality, and the hiPSC‐derived cone photoreceptors are closer to fetal origin compared to their adult counterpart.
In another benchmarking example, we set out to assess the potential differences between
in vitro cell lines compared to adult cells
in vivo. In this regard, we compared the spontaneously immortalized human Müller glia cell line MIO‐M1 (Limb
et al,
2002; Lawrence
et al,
2007) to all the retinal cell types identified in our dataset. Using scRNA‐seq, we profiled 7,150 MIO‐M1 cells with 23,987 reads per cell post‐normalization corresponding to 3,421 detected genes. Unsupervised clustering and t‐SNE analysis showed that the MIO‐M1 cells formed one cluster that is transcriptionally distinct from all retina cell types identified in the human neural retina dataset (Fig
6C). Correlation analysis showed that MIO‐M1 displayed similarities to retinal glial cells, with higher similarity to astrocytes compared to Müller glia (0.63 and 0.46, respectively; Fig
6D). In particular, we identified that MIO‐M1 cells express high levels of the thymosin beta 4 gene (
TMSB4X), which has been linked to glioma malignancy (Wirsching
et al,
2013)
, and the calcyclin gene (
S100A6), which is implicated in macular or cone‐associated diseases (Yoshida
et al,
2004; Fig
5E). Together, our results highlight the similarities and differences in MIO‐M1 to adult retinal glial cells in humans.