Multi-species single-cell transcriptomic analysis of ocular compartment regulons

Identification of cell types in the human and porcine eye

Cells of different regions of the eye were extracted for scRNA-seq, including iris, cornea, choroid, sclera, retina and RPE (Fig. 1a, Supplementary Fig. 1a). In total, approximately 50,000 cells and 24,075 genes were detected. At least 16 distinct clusters were formed by t-distributed stochastic neighbour embedding (tSNE), an unsupervised graph clustering method (Fig. 1b). Annotation of cell types was based on the literature^7,8,9,10 and differential gene profile of each cluster (Fig. 1c and Supplementary Data 1–2). Detailed markers used for cell types are listed in Supplementary Table 1. Characteristic cell type markers were visualised over the tSNE plot to show specific gene expression (Fig. 1d). To add another layer of specificity, we performed a Gene Ontology (GO) analysis of each cluster (Fig. 1e). For photoreceptors (PRs) cells, terms like the sensory perception of the light stimulus were observed, consistent with PRs light-responsive nature. Axon and protein localisation to synapse term was observed in RGCs, accordant with its role in central nervous system connection. Evaluating cell cycle phase genes revealed PRs to have a high proportion of cells in the G2M phase (Supplementary Fig. 1b, e), which might have a role in disc shedding as PRs are known to renew their outer segment by this process¹¹.

One of the advantages of scRNA-seq analysis is the ability to describe the heterogeneity of a given tissue. Here, we observed the RPE layer is relatively homogenous while other tissues have heterogeneous cell populations (Supplementary Fig. 1d). The number of genes expressed per cell is a good indicator of data quality. Checking the quality of cells throughout the tSNE plot by the number of genes expressed indicated overall good cell quality (Supplementary Fig. 1c). As a quality control, we checked for the expression of gender-specific genes and found that they matched the information from donor data (Supplementary Fig. 1f). To probe the quality of our dataset even further, we checked the S-cones and L/M-cones markers in our cone PRs. Both can be detected (Supplementary Fig. 1g) with L/M-cones comprising a much more significant proportion of the whole cone population (Supplementary Fig. 1h). We also checked the phenotype anomalies that might be caused by the mutation of genes in cell types. Phenotype Ontology of rod PR showed retinal dystrophy, while abnormal iris pigmentation was present for melanocytes (Supplementary Fig. 1i). Such analyses will help genetic studies whose mutations contribute to ocular malfunctions. Samples from six human donors were used to create an ocular atlas (Supplementary Fig. 1j). We performed a comparative transcriptome analysis across species and produced ocular atlases for the pig eye (Supplementary Fig. 1k). In total, based on established criteria for scRNA-seq data, the sequence data used here was of good quality.

Diversity and conservation of neural retina across species

Although neural retina on a single cell transcriptomic level has been one of the most well characterised ocular compartments to date (Fig. 2a), new approaches reveal observations not previously understood. In our scRNA-seq data, nine distinct cell types were recognised and sub-classes of bipolar cells (BPs) have been identified (Fig. 2b) based on the canonical markers used for annotation (Supplementary Fig. 2a). First, we checked the similarity of retinal scRNA-seq data with two other published retinal scRNA-seq datasets^5,12. After that, we found a high correlation between corresponding cell types in the eye (Supplementary Fig. 2c), validating our cell classifications. Next, we reported markers with unknown roles in retinal cell types (Supplementary Fig. 2b).

As a validation, RNA fluorescence in situ hybridisation (FISH) was performed using human retinal slides on the candidates of novel markers. Chosen markers NIF3L1 (Fig. 2d) and TRDN (Fig. 2e) were confirmed for localisation in human retinal slides. TRDN was shown to colocalise in VIM + cells in Muller glial cells while NIF3L1 colocalises with PKCA + BPs in RNA FISH experiments.

MALAT1 localisation could be seen in the inner nuclear layer (INL) layer in both targets as an internal control. RAB41 was recently shown to localise in cone PRs¹³. As an additional validation for our RNA FISH results, we offer similar results where RAB41 localises with OPN1SW + cone cells in Fig. 2c. All these novel markers have consistent expression across all the four human donors and may have yet to be explored in terms of functional significance (Supplementary Fig. 2b).

To further explore the NR transcriptome from an evolutionary point of view, we compared our retinal atlas to mouse, primates, and zebrafish retinal atlases^7,14,15 Combinatorial analysis of these datasets enabled the comparison of individual retinal cells across multiple species with the human dataset. We used a single-cell coordinated gene association in pattern sets (scCoGAPS) algorithm to find gene patterns specific to retinal cell types (Fig. 2f). The patterns specific to human retinal cell types were projected into the retina of other species (Supplementary Fig. 2d).

Some pattern-cell type combinations are conserved across species, for instance, pattern 13 for amacrine cells (ACs) (Supplementary Fig. 2e), pattern 24 for cone cells (Supplementary Fig. 2f), pattern 71 for rod cells (Supplementary Fig. 2j), and pattern 6 for rod bipolar (Supplementary Fig. 2k) cells. However, another pattern-cell type like pattern 34 for cone BPs was only found in mammalian species (Supplementary Fig. 2g). Some other patterns of interest were Pattern 4 for horizontal cells (Supplementary Fig. 2h) and Pattern 2 for RGC cells (Supplementary Fig. 2i). All patterns specific to cell types are also listed in Supplementary Data 3. Muller glial cells are of particular interest among all retinal cell types due to their known species difference. The excellent regeneration ability to various neurons in zebrafish marks its significant difference from the mammalian species¹⁶. To understand this, we found patterns specific to Muller glial cells across the species. And among such patterns was Pattern 17 (Fig. 2g), and the genes that constitute Pattern 17 were shown in Fig. 2h. Some of the genes which were included in Pattern 17 were removed (RHO, ACTB, and GAPDH) as they were considered artefacts of the scGoPAS algorithm. We checked the correlation of Pattern 17 to zebrafish Muller glial cells at different stages of activation after injury. We found out that Pattern 17 matches the most with resting Muller glial cells of zebrafish. The similarity decreases as they transit to form progenitors from resting Muller glial cells (Fig. 2i, j). The gene expression patterns that appear in Muller glial cells activation after injury show GO terms like “regeneration”, “sensory system development” and “cell cycle” (Fig. 2k). One similarity of such patterns between zebrafish and humans (like pattern 78 in zebrafish and pattern 17 in humans) was none, as they have very few genes in common (Fig. 2l). Module 17, which was conserved across species, does not include gene patterns related to regeneration.

Characterisation of non-retinal ocular structures

Non-retinal ocular structures help support the structure of the eye and manage light hitting the retina. The sclera provides protection and structure to the eye, while the choriocapillaris (Choroid) is the vascular bed underlying the Bruch’s membrane providing nutritional support for the retina. The single-cell map of the choroidal and sclera layer (Fig. 3a, Supplementary Fig. 3a) shows that cells that populate the sclera are fibrotic tissue. The choroidal layer of the eye consisted of endothelial cells and fibroblast cells. Since it is difficult to separate these two tissues physically and share cell types in common, we analysed them together (Supplementary Fig. 3a). We observed cell types that typically circulate in the choroidal vessels, including activated T cells and monocytes (Fig. 3a). The canonical markers used for annotation of cells and novel markers for those cell types provide additional resources for further study into diseases of the sclera and choroid (Supplementary Fig. 3d).

The cornea is the outermost transparent layer of the eye, whose primary function is transparency and to act as a barrier. The adult cornea has three layers: an outer epithelium (ectoderm), a middle layer containing a collagen-rich stromal region composed of fibroblast cells, and an endothelial cells’ inner layer. Corneal fibroblasts are originated from neural crest cell¹⁶. In our dataset, the cornea is populated mainly by two epithelial cells (Fig. 3b, Supplementary Fig. 3b). One of them have high expression of the TGFBI gene than the rest of the cells, and another has high expression of the ELF3 gene (Fig. 3b, Supplementary Fig. 3e). ANXA1 is a marker of cells undergoing inflammation¹⁷ is expressed highly in TGFBI-high epithelial cells (Supplementary Fig. 3e), while ELF3 is expressed in differentiating corneal epithelial cells¹⁸. Corneal wound healing involves inflammation, proliferation and differentiation processes¹⁹ and expression of TGFBI and ELF3 can distinguish which stage any given corneal cell is present. Such processes might be visualised by the RNA velocity analyses of corneal cells (Supplementary Fig. 3g). Other corneal cell types include fibroblasts, melanocytes, monocytes, cytotoxic T cells, and conjunctival cells. The canonical markers used for annotation of these cell types and the novel markers associated with them are in Supplementary Fig. 3e.

The iris functions to regulate the amount of light reaching the retina. The anterior portion of the optic cup of the eye during development gives rise to iris epithelium and ciliary body epithelium⁸. The stromal region of the iris is generated from neural crest cell migration²⁰. However, the smooth muscles of the iris are developed from optic cup neuroectoderm, and ciliary muscles responsible for changing the shape of the lens are made from surrounding mesenchymal cells⁹. As a result, iris tissue has heterogeneous cells derived from different developmental origins. We found that most cells were fibroblasts from our scRNA-seq data of the iris/ciliary body (Fig. 3c, Supplementary Fig. 3c). They were subdivided into fibroblasts, MEG3-high fibroblasts, MGP-high fibroblasts, WIF1-high fibroblasts and Ribosomal genes high fibroblasts. Fibroblasts highly expressing MEG3 is proliferating in glaucoma tenon fibroblast¹⁰.

Similarly, fibroblasts having high expression of WIF have been shown to initiate melanogenesis in normal human melanocytes²¹. Iris stroma is also populated by Schwann cells that help myelination of axons of neuronal cell types that populate iris²². They expressed CD9 marker, which is implicated in the signalling of Schwann cells with axon²³.

Ciliary body cells (CBCs) help produce aqueous humour in the anterior chamber of the eye²⁴. CBCs were subdivided into CBCs, COL9A1- high CBCs, CRYAA-high CBCs and pigmented CBCs (Fig. 3c). Smooth muscle cells (SMCs) that populate sphincter pupillae and dilator pupillae muscles were detected in the data (Fig. 3c). They had high MYH11 and MYL9 (Supplementary Fig. 3f) expression, which are the canonical markers for SMCs²⁵. A pigmented layer of epithelial cells also populates the iris. They expressed canonical markers involved in  pigmentation like MLANA and PMEL (Supplementary Fig. 3f).

We tried to visualise the cell-cell interaction among all cell types of the eye. As expected, we found more interactions with cells located physically together (Supplementary Data 4). However, monocytes were also seen to have more interactions with the cell types in the retina. We found the most significant interacting pairs (Fig. 3d) and clustered them using the tSNE plot (Fig. 3e). The tSNE plot helps to give an idea about the similarity of cell–cell interactions among cell types. We aimed to provide a global picture of cell–cell interactions across cell types, and Fig. 3f showed these kinds of interactions.

Also, we checked cell types having the most interactions in the whole eye. The most interactions were with Muller glial cells across species (Supplementary Fig. 3h). However, such analysis is limited by the genes selected for the study, which were paralogs of human genes across species. As a result, we could see lesser interactions in other species compared to humans.

Putative stem cell populations of the iris

Among various cell types in the iris/ciliary body, we could detect a distinct population of cells that express markers of stem cells. Even though these cells express markers of retinal progenitor cells or retinal stem cells, such cells need to give rise to retinal neurons in vivo to label them as retinal stem cells. The absence of in vivo experiments limits our work. Previous studies have mentioned that Muller glial cells share 60% of the transcriptome with retinal progenitors in mice²⁶. However, transcriptomic similarity or the presence of specific genes markers do not guarantee if a cell is a retinal stem cell or not. Therefore, such cells remain as putative stem cells for now.

Besides putative stem cells, two subpopulations of CBCs, COL9A1 high CBCs and pigmented CBCs showed higher similarity to such putative stem cells. All three cell types expressed of PAX6 and SIX3, eye field TFs (EFTF). Putative stem cells have an expression of OTX2, which is another EFTFs (Supplementary Fig. 3f). COL9A1 high CBCs expressed CPAMD8 (Supplementary Fig. 3f), which plays a role in periocular mesenchyme development²⁷. We also found the presence of both cells in the iris of pigs (Fig. 4c). Checking for stem cell potency across cell types showed the identified stem cell populations score higher on a stem-cell potency index (Fig. 4a, b), providing evidence that such cell populations may exhibit stem cell properties. RNA velocity is a high dimensional vector that predicts the fate of cell populations in scRNA-seq data²⁸. Here, RNA velocity analysis showed all three putative stem cells cluster together in Uniform Manifold Approximation and Projection (UMAP) plots (Supplementary Fig. 4a). It is suggestive of transcriptional similarity hinting that pigmented CBCs and COL9A1-high CBCs might originate from putative stem cells. We tried to understand the gene patterns that make up the putative stem cells, pigmented CBCs and COL9A1-high CBCs. We found out that patterns 13, 37, 48, 54, 55, 56, 57, 86 and 87 was shared among three cell types (Fig. 4d). When these gene patterns were projected into pig iris cell types (Supplementary Fig. 4b), we could see that specific gene patterns like patterns 13, 87, 56 and 86 were conserved across the species (Fig. 4d, e). Upon further examining the type of biological processes these gene patterns constitute, we investigated the GO terms that were enriched in such gene patterns. GO terms like “Neural Crest Differentiation”, Negative regulation of cell differentiation”, and “embryonic morphogenesis” were enriched (Fig. 4f). One of the patterns specific to COL9A1 high CBCs, Pattern 56, was probed for the genes it comprised (Fig. 4g). Previous research has pointed that there are multipotent cells derived from neural crest in adult mouse iris²⁹. So, this suggests that COL9A1 high CBCs might be some sort of multipotent stem cells.

Ligand–receptor interactions in putative ocular stem cells

We focused our analysis on signalling pathways involved in stem cell maintenance to identify which primary molecules are at play in the eye. We focused on three signalling pathways, i.e., Fibroblast growth factor (FGF), WNT and Midkine (MDK) signalling pathways (Fig. 4i–k). Midkine, one of the ligands specific to retinal progenitor cells in zebrafish, showed high expression in putative stem cells³⁰. It was also shown to mediate glial activity, neuronal survival and the reprogramming of Muller glia into proliferating Muller glial proliferating cells in chicks³¹. MDK expression was high in putative stem cells, and its receptors were shown to be present in both neuronal and non-neuronal cell types of the eye (Fig. 4k). Thus, the expression of MDK added to the evidence of putative stem cells being present in iris tissue. When we focussed on the WNT signalling pathway, gene expression of ROR1, ROR2 and FZD1 receptors were found to be distributed among putative stem cells, pigmented CBCs and COL9A1 high CBCs (Fig. 4h). MAGIC co-expression analysis based on MAGIC imputation revealed WNT receptors ROR1, ROR2, and FZD1 were specific to pigmented CBCs, putative stem cells and COL9A1-high CBCs, respectively. These cells were also uniquely co-expressed with FGFR1 (Supplementary Fig. 4e). FGFR1 expression was also conserved across similar cell types in pigs (Supplementary Fig. 4c, d).

Creation of disease map and viral-entry map for the human eye

We wanted to provide a resource to understand the disease map and viral-entry map of cell types in the human eye. For that, we obtained a list of genes that were affected in ocular malformations and checked the gene expression of those genes across the whole eye. We focused on genes that cause colour blindness, corneal disorders, eye cancer, eye movement disorders, macular degeneration, optic nerve disorders, retinal disorders, vision impairment and blindness (Fig. 5a). For example, GSN mutation causes reduced corneal sensitivity in the later stage of life. Mutations in GSN causes deposition of amyloid in a different part of the eye³² However, expression of GSN throughout the cell types of the eye shows that it is expressed in CBCs, SMCs and fibroblasts in the eye (Fig. 5b). It gave an idea of how mutations of GSN could cause corneal dystrophy.

We also checked the expression of such genes across species. We showed that genes for cone/rod dystrophy, retinitis pigmentosa and stationary night blindness are conserved across species (Supplementary Fig. 5a). We looked into the viral entry map of the human eye. As the world’s interest in COVID-19 increases, we wanted to provide a map of cell surface proteins that can act as viral entry receptors. We showed that ACE2 and TMPRSS2, the primary cell surface proteins responsible for entry into human³³, are expressed in the cornea’s conjunctival cells (Supplementary Fig. 5b, c). Other receptors of interest like BSG, CTSB and CTSL were also shown. As a resource, we provided information for the cell surface proteins that serve as an entry point for other classes of viruses (Supplementary Fig. 5b).

Unique transcriptional regulons active in the human eye

Since TFs orchestrate gene expression across the genome, cell identity and function can be partially described by the expression of its TFs. The TFs expressed may provide insight into the machinery that maintains their stemness with a focus on the putative stem cells. We made a pipeline to understand transcriptional regulons active in different retinal and putative stem cells. We combined the Regulon activity score computed from SCENIC and gene imputation scores calculated from MAGIC to create a pairwise correlation of TFs active in cell types (Supplementary Fig. 6a). The TFs which were enriched in cell types are also listed in Supplementary Data 5. We found 9 modules specific to cell types in the retina (Fig. 6a, Supplementary Fig. 6b). Module M1 was specific to RGCs, and module M5 was specific to PR cells and presented in Fig. 6b, c, which included the responsible TFs. The specificity of such modules in cell types could be seen in tSNE plots (Fig. 6d). With the help of a correlation matrix among TFs, we plotted the interactions of TFs with each other in modules (Fig. 6e).

Analysis of the correlation map showed that TFs active in Schwann cells and Muller glial cells were highly correlated as both are glial cell types. The pigmented CBC types and retinal pigmented epithelial cells were also highly correlated as both are pigmented cells. Besides checking the transcriptional modules of retinal cells, we evaluated the transcriptional modules of non-retinal cell types (Supplementary Fig. 6c–e) of the eye. The TFs that populate melanocytes of Iris included PAX3, MITF and SOX10 (Supplementary Fig. 6e). These factors have been shown to be important in the trans-differentiation of fibroblasts into melanocyte³⁴. Thus, such TF modules might also provide a resource for trans-differentiation of cell types.

Conservation of TF modules across species

We checked for the genes that were involved in the formation of TF modules. To do so, we conducted a GO analysis of the TFs, and their target genes confirmed the GO terms are specific to each cell type (Fig. 7a, Supplementary Data 6). To gain additional confidence in the presence of TF module specific to cell types, we verified SREBP2 and KLF7 specificity in RGCs and PBX1 in ACs and RGCs using RNA FISH. SREBP2 and KLF7 are two TFs present in Module 1 specific to RGCs (Fig. 6e).

SREBP2 is TF involved in cholesterol biosynthesis³⁵ and KLF7 involved in the axon regeneration response after optical nerve injury in the eye³⁶. We detected the localisation of SREBP2 and KLF7 in RGC cells (Fig. 7b, Supplementary Fig. 7a). PBX1 localisation to the INL layer where amacrine cells were located could also be observed. RLBP1 was used as a control where it localised to the RPE layer (Supplementary Fig. 7a). Besides RNA FISH, we also conducted immunofluorescence in the primate retina for the TF specificity. We checked whether genes involved in TF modules are sufficient to separate cell types across species and PCA analysis and demonstrated clear separations between cell types based on the unique TFs (Supplementary Fig. 7b, c) compared to randomly selected some TFs (Supplementary Fig. 7d, e). Such conservation of TF modules could be shown across species (Supplementary Fig. 7f) using pairwise correlation. The TFs that were involved in the formation of modules were checked for their conservation. For example, pairwise correlation plots show that LHX9, TFAP2C were conserved across species for amacrine cells.

Similarly, KLF7 was shown to be preserved in RGCs across humans, mice, macaque and zebrafish (Supplementary Fig. 7f). Such analysis helped us to understand the TF differences in cell types across the species. In Klf7-null mice, a small portion of RGCs showed aberrant projections while exiting retina³⁷. Klf7 and Pou4f1, another member of RGC module TF, co-operate to control TrkA expression in sensory neurons³⁸. Since KLF7 is also expressed in foetal RGC during RGC development in humans³⁹, conserved across species and has a role in axon regeneration response after injury, we wanted to focus on the part of KLF7 in RGC differentiation and maturation.

KLF7 acts as a driver for RGC maturation

We used a protocol⁴⁰ to differentiate H9 Human embryonic stem cells into RGCs. We performed KD experiments in cells that were driven to RGC lineage (Fig. 7c, d, Supplementary Fig. 7g, h). Successful KD was achieved by shRNA transfection (Fig. 7e). Since cells were still in the process of differentiation, mature RGC markers were lowly expressed. However, markers like POU4F1 and EBF3, which act as markers for maturing RGC^41,42, were present in cells (Supplementary Fig. 7j). Cells expressing ONECUT1, ONECUT2 and EBF3 were designated as RGC precursor cells because ONECUT TFs are expressed in developing RGCs⁴³.

As a result of the KD of KLF7, we detected drastic differences in the proportion of cells that are destined to be RGCs (Fig. 7f). Also, genes like POU4F1, IRX2 and EBF3, which comprised the TF modules of RGC (M1 module), were downregulated after KD. Moreover, in RGC precursor cells, decreased expression of EBF1, an M1 module TF, was observed (Fig. 7j). To gain further insights, we overexpressed KLF7 during the early differentiation window (Fig. 7g, Supplementary Fig. 7i), which aimed to verify whether the acceleration of RGC differentiation can be achieved by KLF7 alone. OE was confirmed (Fig. 7h) after manual annotation with the same criteria used for KD libraries (Supplementary Fig. 7k). Expectedly, the early differentiation window did not yield many mature RGCs. However, the proportion of RGC-like cells was increased after KLF7 OE (Fig. 7i). Moreover, KLF7 OE also increased the expression of RGC module TFs, like EBF3, SREBF2 and EBF1 in RGC precursor cells (Fig. 7k). Those changes hint that KLF7 might have a role in biasing cell fate toward RGCs during retinogenesis.