Single cell transcriptomics reveals the heterogeneity of the human cornea to identify novel markers of the limbus and stroma

The cornea is the clear window that lets light into the eye. It is composed of five layers: epithelium, Bowman’s layer, stroma, Descemet’s membrane and endothelium. The maintenance of its structure and transparency are determined by the functions of the different cell types populating each layer. Attempts to regenerate corneal tissue and understand disease conditions requires knowledge of how cell profiles vary across this heterogeneous tissue. We performed a single cell transcriptomic profiling of 19,472 cells isolated from eight healthy donor corneas. Our analysis delineates the heterogeneity of the corneal layers by identifying cell populations and revealing cell states that contribute in preserving corneal homeostasis. We identified expression of CAV1, HOMER3 and CPVL in the corneal epithelial limbal stem cell niche, CKS2, STMN1 and UBE2C were exclusively expressed in highly proliferative transit amplifying cells, CXCL14 was expressed exclusively in the suprabasal/superficial limbus, and NNMT was exclusively expressed by stromal keratocytes. Overall, this research provides a basis to improve current primary cell expansion protocols, for future profiling of corneal disease states, to help guide pluripotent stem cells into different corneal lineages, and to understand how engineered substrates affect corneal cells to improve regenerative therapies.

The healthy cornea is a transparent and avascular tissue that allows light to enter the eye and accounts for most of its refractive power. The cornea is composed of five layers: its outer surface is a stratified sheet of corneal epithelial cells that reside on the Bowman's layer, a collagen-based acellular membrane synthesized by the stromal keratocytes. The keratocytes populate the corneal stroma, the middle layer of the cornea that is composed of structured collagen fibers and other extracellular matrix proteins. The corneal endothelium forms a thin monolayer of tightly packed hexagonal cells that line in the innermost surface of the cornea and reside in close contact to the stroma on the Descemet's membrane.
Corneal structure and transparency are governed by the functions of the cell types populating each layer. Epithelial cells act as a biological barrier to block the passage of foreign material and provide a smooth surface that absorbs nutrients. Keratocytes maintain extracellular matrix homeostasis responsible for the cornea's biomechanical and optical properties, and endothelial cells serve as active pumps transporting ions, metabolites and fluid to maintain corneal hydration and transparency.
Selective replacement of dysfunctional single corneal layers with that of a donor 1 , the autologous transplantation of primary cultured corneal epithelial limbal stem cells 2 , and the treatment of endothelial dysfunctions with primary cultured allogeneic corneal endothelial cells 3 are already a therapeutic reality. Gaining deeper understanding of corneal cell profiles and their transcriptomic signatures can be highly relevant for the improvement of such therapies.
In order to further understand the cellular complexity of this heterogeneous tissue, we provide a high-quality single-cell ribonucleic acid sequencing (scRNAseq) dataset from 19,472 corneal cells isolated from four female and four male donors. With it, we depict the diverse cell populations and provide a comprehensive cell atlas of Nine cell clusters were identified within the corneal epithelium. A major cluster of 5964 corneal epithelial cells was identified and further analysis revealed nine subclusters (E0-8, Fig. 3a). Differential gene expression profiling was used for further identification of the subclusters. Cluster E3 presented a high expression of stromal cell markers KERA, lumican (LUM), and aldehyde dehydrogenase 3 member A1 (ALDH3A1) 5,7,8 , as well as a reduced expression of corneal epithelial marker keratin 12 (KRT12) 9,10 compared to the other epithelial subclusters ( Fig. 3b and Supplementary Figure S1). Putative doublet analysis of cluster E3 suggested epithelialstromal keratocyte cell doublets (Supplementary Figure S2). Cluster E7 was identified as conjunctival epithelial cells based on high expression of conjunctival markers such as keratin 13 (KRT13), keratin 15 (KRT15), and keratin 19 (KRT19) 11,12 , and the low expression of the corneal epithelial marker KRT12 ( Fig. 3b and Supplementary Figure S1).
Cells forming clusters E6, E8 and E8 showed increased expression of corneal limbal markers keratin 14 (KRT14) 13,14 , KRT15 15,16 , and S100 calcium binding protein A2 (S100A2) 17 compared to the other identified clusters suggesting their location in the corneal limbus ( Fig. 3b and Supplementary Figure S1). Interestingly, cluster E6 showed an increased expression of the superficial epithelium limbus marker 17 S100 calcium binding protein A8 (S100A8) compared with clusters E0 and E8 and a reduced expression of the basal corneal epithelial cell markers connexin 26 (GJB2), connexin 30 (GJB6), and integrin β4 (ITGB4) 18-20 , (Fig. 3b and Supplementary Figure S1) leading to its identification as a population of wing/superficial epithelial cells in the limbus or peripheral cornea. Cluster E0 presented an increased expression of GJB6, and GJB2, which are predominantly found in basal corneal epithelium 18 suggesting these cells formed a basal corneal epithelial cell population in the limbal stem cell niche or peripheral cornea ( Fig. 3b and Supplementary Figure S1). Finally, the high expression of mitogenic factors Ki-67 (MKI67) 21 , survivin (BIRC5) 22 and H2A histone family member X (H2AFX) 23 in cluster E8, as well as the differential expression of transit-amplifying cell marker CD109 24 , suggested that cluster E8 was formed by highly proliferative transit-amplifying cells in the limbal stem cell niche or peripheral cornea ( Fig. 3b and Supplementary Figure S1). No quiescent limbal epithelial stem cells expressing ABCB5, ABCG2 and CD200 [24][25][26][27] were identified in this dataset (Supplementary Figure S1).
Cluster E5 presented high expression of GJB6 and ITGB4 associated with basal epithelium 18 and KRT12, and KRT3 associated with terminally differentiated corneal epithelium [27][28][29] , along with a reduced expression of KRT14  Figure S1). Cluster E1 was identified as post-mitotic and terminally differentiated migratory epithelial cells based on the expression of genes associated with cell migration such as RHOV 30 and tight junction formation and obliteration CLDN7 31 together with a high expression of corneal epithelial cell markers KRT12, KRT3, and KRT5. This cluster retained a low expression level of KRT14, implying their limbal origin ( Fig. 3b and Supplementary Figure S1). Cluster E4 was composed of cells with a high expression of KRT12, KRT3, and KLF5 suggesting this cluster was terminally differentiated cells from the central corneal epithelium (Fig. 3b and Supplementary Figure S1). The cells forming cluster E2 presented a high expression of KRT12, KRT24, and CXCL17 associated with wing/ superficial central epithelium ( Fig. 3b and Supplementary Figure S1). scRNAseq reveals novel specific markers for the corneal limbal stem cell niche and transit-amplifying cells. Differential gene expression analysis of the basal corneal limbal epithelial cells (cluster E0), and the transit-amplifying cells (cluster E8), showed that genes encoding for caveolin-1 (CAV1), probable serine carboxypeptidase (CPVL), homer scaffolding protein 3 (HOMER3), and C-X-C motif chemokine 14 (CXCL14) were highly expressed in both clusters (Fig. 3c), suggesting they could be markers of the human limbal stem cell niche. Moreover, the expression of the markers highly correlated in clusters E0 and E8 (Fig. 3d), again To confirm these findings, expression of the markers in the corneal epithelial limbus and central cornea was assessed by immunofluorescence ( Fig. 4 and Supplementary Figure S3). Caveolin-1 and CXCL14 were expressed in the limbus of the cornea, whereas caveolin-1 was minimally expressed and CXCL14 absent in the central cornea ( Fig. 4 and Supplementary Figure S3 Figure S4).
Differential gene expression analysis on the transit-amplifying cells (cluster E8) revealed that the expression of CKS2, STMN1, and UBE2C was exclusive to this cluster. Expression of cyclin-dependent kinase 2, stathmin-1, and ubiquitin conjugating enzyme E2 C, was exclusive to the limbus and absent in the central cornea (   Fig. 5a). High expression of stromal keratocyte markers LUM, KERA, or ALDH3A1 5,7,8 suggested that clusters S0, S1 and S2 were stromal keratocytes (Fig. 5b). Cluster S2 showed increased expression of extracellular matrix proteins such as COL12A1, COL6A2, COL6A1, and LAMB2 suggesting these were activated keratocytes that played a crucial role in maintaining the corneal stromal extracellular matrix ( Fig. 5b and Supplementary Figure S7). Cluster S3 showed a decreased expression of keratocyte markers LUM, KERA, and ALDH3A1 and increased expression of the fibroblastic marker CD44 32-34 compared to clusters S0, S1 and S2 ( Fig. 5b and Supplementary Figure S7). The myofibroblast-specific marker, α-smooth muscle actin (ACTA2) 35,36 , was not detected in cluster S3 (Supplementary Figure S7), suggesting it is composed of keratocytes that are in transition to stromal myofibroblasts. Finally, differential gene expression analysis of the major stromal cluster compared to the clusters associated with the corneal epithelial and endothelial layers identified that the expression of nicotinamide N-methyltransferase (NNMT) was exclusive to the corneal stromal cluster, suggesting that this gene could be used as a novel corneal stromal marker (Fig. 5b,c).  Fig. 6a). Both clusters showed expression of corneal endothelial cell markers CD166 (ALCAM) and sPrdx1 (PRDX1) 37,38 and functional markers Na + /K + ATPase (ATP1A1) 6 and sodium bicarbonate transporter-like protein 11 (SLC4A11) 6 confirming the corneal endothelial phenotype ( Fig. 6b and Supplementary Figure S7). Differential expression analysis revealed that cluster En0 possessed a lower expression of tight junction protein zona occludens-1 (TJP1) 6 and focal adhesion regulator microtubule-actin cross-linking factor-1 (MACF1) 39 compared to cluster En1 (Fig. 6b), suggesting cells in cluster En0 could preferentially migrate upon corneal endothelial damage to contribute to tissue repair. Furthermore, cluster En1 possessed a higher expression of COL4A3 (Fig. 6b) suggesting these cells could play an important role on Descemet's membrane homeostasis. Interestingly, both clusters retained expression of PITX2, a periocular mesenchyme marker associated with progenitor endothelial cells 40,41 . No endothelial fibroblasts were identified, evidenced by the absence of ACTA2 and CD44 (Supplementary Figure S8).

Discussion
Gaining transcriptomic information at the single cell level of human corneal cells enables a greater understanding of this heterogeneous tissue. In this study, we have performed a scRNAseq analysis of the healthy cornea to create a comprehensive cell atlas of the human cornea (Fig. 7). Moreover, scRNAseq analysis enabled the identification of novel markers of the limbal epithelial stem cell niche, transit-amplifying cells, and stromal keratocytes. The data generated can serve as a reference cell atlas with a major impact in the further improvement and development of cell replacement therapies or regenerative medicine approaches for treating corneal The corneal epithelium appeared to be the most heterogeneous corneal layer with 9 identified cell clusters. The transcriptomic signature of epithelial cells from the basal, wing and superficial layers of both central and limbus/peripheral cornea were identified. A population of conjunctival epithelial cells was detected (cluster E7; Fig. 3), likely representing contamination from the dissection process. Moreover, a cluster of highly proliferative  Differential analysis identified the expression of CAV1, HOMER3, CXCL14, and CPVL to be exclusive to cell clusters comprising the corneal epithelial stem cell niche (Fig. 3). In line with our findings, Collin et al. recently reported the exclusivity of CXCL14 and CPVL to the corneal epithelial stem cell niche 42 . Immunofluorescence analysis confirmed the expression of caveolin-1 and CXCL14 in the limbus, and their absence in the central cornea (Fig. 4). Nevertheless, CXCL14 was expressed in the superficial layers of the limbus but not the limbal stem cell niche. HOMER3 and CPVL were both expressed in the limbus, and in the central cornea, where central basal epithelial cells appeared to retain higher expression of these markers (Fig. 4). We further validated the finding of caveolin-1 with immunofluorescence in primary cultured limbal stem cells (Supplementary Figure S4). These results suggest that caveolin-1 could be used for the identification of epithelial limbal stem cells, with the , and a cluster of highly proliferative transit-amplifying cell (E8). In the stroma: three keratocyte clusters (S0-S2), with one having a high extracellular matrix protein secretion profile (S2) and a cluster of keratocytes transitioning to stromal myofibroblasts (S3). In the endothelium: a cell cluster with high collagen synthesis (En1), and a cell cluster with low tight junction and focal adhesion protein expression (En0). www.nature.com/scientificreports/ advantage that it is a cell membrane marker, whereas p63 is expressed in the nucleus. This could open the door to isolating and enriching these cells for regenerative therapies. Furthermore, our study also revealed novel makers specific to transit-amplifying cells, namely CKS2, STMN1, and UBE2C. We further validated these findings with immunofluorescence, and confirmed their expression in the limbus and periphery of human corneas. Furthermore, expression of CKS2 was also assessed in primary cultured limbal stem cells (Supplementary Figure S6), suggesting this is a suitable marker for identifying highly proliferative transit-amplifying cells. Interestingly, primary limbal cells expressing a lower level amount of ΔNp63 and p63α, based on fluorescence intensity, also showed a reduced expression of caveolin-1 and cyclin-dependent kinase 2 respectively (Supplementary Figures S4  and S6). Interestingly, two recent studies also reported STMN1 and UBE2C to be specific markers of limbal transit amplifying cells 43,45 . Finally, in a 2019 single cell transcriptome study, Kaplan et al. hypothesized that the stemness of mouse limbal epithelial stem cells could be regulated through autophagy 46 . The identification of caveolin-1, a facilitator of caveolin-mediated endocytosis, as a marker of the human corneal limbus in our study could support this hypothesis. Nevertheless, further research is required to confirm this theory.
It is important to highlight that as previously reported, the corneal epithelium sheds superficial layers when corneas are preserved in media 47 and it is very likely that this cell atlas does not fully portray the most superficial layer of corneal epithelial cells. Two types of cornea preservation conditions, namely in Optisol-GS medium and organ culture medium, were used for this research, in line with clinical practice in the United States of America and Europe, respectively. No major differences in cell clustering were found between preservation conditions. Three clusters of the corneal stroma were identified as corneal keratocytes (clusters S0, S1 and S2; Fig. 5). Cluster S2 was identified as cells with a key role in extracellular protein secretion to maintain the homeostasis of the stroma. Furthermore, we identified a cluster (S3) of keratocytes transitioning to myofibroblasts, though a fully differentiated myofibroblast phenotype expressing alpha smooth muscle actin was not identified. Interestingly, our analysis did not identify a population of corneal stromal stem cells expressing ABCB5 and ABCG2, as previously identified by Funderburgh and colleagues [48][49][50] . Our hypothesis is that the expression of these genes is induced by the primary expansion of corneal keratocytes and is not found in the cornea [48][49][50] . Finally, differential gene expression analysis identified that the expression of NNMT was exclusive to the corneal stroma, suggesting this could be used as a novel marker to identify corneal stromal cells. Interestingly, the study by Li et al. reported differential expression of NNMT in a cluster of 69 quiescent limbal stem cells expressing ABCB5 and ABCG2 43 . As previously mentioned, perhaps due to limitations in sequencing depth, we did not find a population of limbal stem cells expressing ABCB5 and ABCG2, and therefore unable to confirm these results. Since Li and colleagues did not include a population of corneal stromal cells in their differential analysis, it also explains why NNMT was not detected as a stromal marker in their study. It is possible that NNMT is expressed specifically in both ABCB5 + /ABCG2 + quiescent limbal stem cells and in corneal stromal cells.
Two cell clusters were identified forming the corneal endothelium (Fig. 6). Both clusters showed high expression of corneal endothelial markers such as CD166 (ALCAM), sPrdx1 (PRDX1), ZO-1 (TJP1), or SLC4A11, confirming the endothelial phenotype. Nevertheless, expression of CD200, reported as corneal endothelial cell marker in a previous study 51 , was not detected in our dataset, suggesting it might not be a specific marker for corneal endothelial cells. Despite sharing great similarity, one of the clusters (En1) expressed more extra cellular matrix proteins, suggesting that these cells could play a crucial role on maintaining the Descemet's membrane. Cell cluster En0 showed lower expression of tight junction proteins and focal adhesions, suggesting these cells could preferentially migrate upon corneal endothelial damage and contribute to tissue repair via migration or cytosolic expansion. Interestingly, our study did not detect a population of precursor endothelial cells, as discussed in previous studies 52,53 . Nevertheless, both corneal endothelial cell clusters detected retained expression of PITX2, a marker associated with neural crest-derived corneal endothelial cell precursors 40,41 . These data suggests that it is highly possible that a progenitor-like state is not exclusive to the peripheral corneal endothelium but to cells across the endothelium.
Interestingly, no corneal endothelial fibroblasts were detected in this study, indicated by the lack of expression of fibroblastic markers α-smooth muscle actin (ACTA2) and CD44 (Supplementary Figure S8), which is in contrast to a previous report 42 . It is likely that the bulk enzymatic corneoscleral tissue desegregation performed by Collin et al., affected the susceptible corneal endothelial cells, causing endothelial cells transcriptomic bias portrayed in scRNAseq dataset, as showed in other studies 54,55 . We performed a more gentle approach to obtain single endothelial cells for sequencing, first dissecting the tissue and mechanically stripping the endothelium from the cornea, and then treating the endothelial cells with a shorter enzymatic digestion.
Overall, this study provides significant information to help understand the heterogeneity of the healthy human cornea, as well as understanding the gene expression stratification across cells present in the same corneal layer, while providing novel markers to identify specific cell types. Moreover, this transcriptomic cell atlas offers a baseline for future studies with the aim of regenerating corneal tissue or further developing corneal cell replacement therapies.

Materials and methods
Ethical statement. This study was performed in compliance with the tenets of the Declaration of Helsinki.
Ten human donor corneas (Table 1)  Manual dissection of the corneal tissue. Eight donor corneas were dissected for scRNAseq. First, cells were isolated from six corneas that were manually dissected to separate the epithelial, stromal and the endothelial layers to ensure the gentlest possible enzymatic dissociation of each layer. The corneas were vacuum fixed in a punch base (e.janach) endothelial-cell side up, stained with trypan blue solution (0.4%) for 30 s and washed with BSS ophthalmic irrigation solution. The corneal endothelium was then gently lifted following the Schwalbe line using a DMEK cleavage hook (e.janach) and fully stripped using angled McPherson tying forceps. The remaining tissue was trephined using a 9.5 mm Ø Barron vacuum punch (Katena) in order to separate the epithelium and stroma from the scleral ring. The remaining two corneas were used for limbus isolation. A surgical scalpel was used to cut the limbus into approximately 1 × 2 mm fragments, which were then rinsed with PBS.
Tissue dissociation to single cells. The six manually dissected corneal tissues were enzymatically treated to obtain single cell suspensions. The stripped corneal endothelium was incubated with 2 mg/mL collagenase (Sigma) solution in human endothelial serum free media (SFM) (Thermo Fisher Scientific) for 1-2 h at 37 °C followed by a 10 min incubation with Accutase (Thermo Fisher Scientific) at 37 °C to obtain a single cell suspension. The cells were centrifuged for 5 min at 800×g and resuspended in 0.5 mL of human endothelial SFM media. The corneal epithelium-stroma tissues were treated with 1.5 mg/mL collagenase and 0.2 mg/mL bovine testis hyaluronidase (Sigma) in DMEM/F12 (Thermo Fisher Scientific) for approximately 1 h at 37 °C. The released corneal epithelium was suctioned with a P1000 micropipette and further treated with Accutase for 10 min at 37 °C to obtain a single cell suspension of corneal epithelial cells. The cells were centrifuged for 5 min at 800×g and resuspended in 0.5 mL of DMEM/F12. The remaining corneal stroma was treated with 1.5 mg/mL collagenase and 0.2 mg/mL bovine testis hyaluronidase in DMEM/F12 for 5-7 h to obtain a single cell suspension of corneal keratocytes. After the incubation, the cells were centrifuged for 5 min at 800×g and resuspended in 0.5 mL of DMEM/F12 media.
The limbus fragments isolated from the two corneas for limbal isolation were digested with four trypsinization cycles. In each cycle, the limbus fragments were incubated in 10 mL of 0.05% trypsin/0.01% EDTA (Thermo Fisher Scientific) at 37 °C for 30 min. The trypsin containing dissociated cells were transferred to a 50 mL centrifuge tube containing 10 mL DMEM/F12 media and centrifuged for 5 min at 300×g, after which the cells were resuspended in 0.5 mL of DMEM/F12 media. The undigested limbal fragments were placed again in 10 mL of 0.05% trypsin/0.01% EDTA for the following trypsinization cycle. The cells isolated from the limbus of the two donor corneas were pooled for the following steps.
Methanol cell fixation. After dissociation, the single cell suspensions were passed through a 100 μm Ø cell strainer, centrifuged for 5 min at 300×g and resuspended in 1 mL ice-cold PBS to eliminate any medium remnants. Next, the cells were again centrifuged for 5 min at 300×g and resuspended in ice-cold PBS at a ratio of 200 μL PBS per 1 × 10 6 cells followed by the dropwise addition of ice-cold methanol, at a ratio of 800 μL per 1 × 10 6 cells. The cells were stored at − 80 °C until sequencing.
Single-cell RNA sequencing (scRNAseq). Single-cell mRNA sequencing was performed at Single Cell Discoveries according to standard 10 × Genomics 3' V3.1 chemistry protocol. Prior to loading the cells on the 10 × Chromium controller, cells were rehydrated in rehydration buffer. Cells were then counted to assess cell integrity and concentration. Approximately 9000 cells (comprising 3000 cells from each of the three layers) from corneas 01, 02, 03 and 04 were separately loaded by layer and cornea. Furthermore, 3000 cells (1000 cells from each of the three layers) from corneas 05 and 06, and 3000 cells from corneas 07 and 08 (from the limbal samples) were separately loaded by cornea. The resulting sequencing libraries were prepared following a standard 10 × Genomics protocol.
Bioinformatic analysis of scRNA-seq data. BCL files resulting from sequencing were transformed to FASTQ files with 10 × Genomics Cell Ranger mkfastq. FASTQ files were mapped with Cell Ranger count. During sequencing, Read 1 was assigned 28 base pairs, and were used for identification of the Illumina library barcode, cell barcode and UMI. R2 was used to map the human reference genome GRCh38. Filtering of empty barcodes was done in Cell Ranger. The data from all samples were loaded in R (version 3.6.2) and processed using the Seurat package (version 3.2.0) 56 . More specifically, cells with at least 1000 UMIs per cell and less than 20% mitochondrial gene content were retained for analysis. The data of all 10 × libraries was merged and processed together. The merged dataset was normalized for sequencing depth per cell and log-transformed using a scaling factor of 10,000. The patient effect was corrected using the integration function of Seurat and used for dimensionality reduction and clustering of all cells or cells selected per layer. Cells were clustered using graph-based clustering and the original Louvain algorithm was utilized for modularity optimization. The differentially expressed genes per cluster were calculated using the Wilcoxon rank sum test and used to identify cell types. Putative doublets were computationally identified using scDblFinder (v1.2.0) 57 .

Primary culture of limbal cells.
Human primary limbal cells were harvested from corneal tissue of cadaveric donors (ages ranging from 36 to 79 years) with informed consent. Human limbal cells were cultured as previously described 58,59 . In short, a surgical scalpel was used to cut the limbus of the corneas into approximately 1 × 2 mm fragments, which were then rinsed with PBS. The limbus fragments were then incubated in 10 mL of www.nature.com/scientificreports/ 0.05% trypsin/0.01% EDTA (Thermo Fisher Scientific) at 37 °C for 30 min. The trypsin containing dissociated cells were transferred to a 50 mL centrifuge tube containing 10 mL culture media and centrifuged for 5 min at 300×g, after which the cells were resuspended in 0.5 mL of culture media. The undigested limbal fragments were placed again in 10 mL of 0.05% trypsin-EDTA for another trypsinization cycle. Culture medium consisted of 2:1 mixture of DMEM/F12 media (Thermo Fisher Scientific) supplemented with 2 mM GlutaMAX (Thermo Fisher Scientific), 10% fetal bovine serum (Thermo Fisher Scientific), 125 IU/L insulin (Humulin R, Lilly), 0.2 mM adenine (Merck), 1.1 μM hydrocortisone (Merck), 8.5 ng/mL cholera toxin (Sigma), 2 nM triiodothyronine (Sigma), 10 ng/mL epidermal growth factor (Amsbio), and 100 IU/mL penicillin-streptomycin (Thermo Fisher Scientific). The limbal cells were plated on a feeder layer of lethally irradiated 3T3-J2 fibroblasts (fibroblast feeder-layer density 40,000 cells/cm 2 ). The 3T3-J2 fibroblast immortalized cells line was a kind gift of Prof. Howard Green (Harvard Medical School, Boston, MA, USA). When confluent, limbal epithelial stem cells were passaged by 0.05% trypsin/0.01% EDTA (Thermo Fisher Scientific) treatment and seeded at a density of 15,000 cells/cm 2 in a Nunc chamber slide (Thermo Fisher Scientific).
Immunofluorescence. Two human donor corneas and primary cultured limbal cells were used for immunofluorescence analysis. The corneas were cut in half transversally, embedded in a cryomold containing Tissue- Tek O.C.T. compound, snap frozen in liquid nitrogen, and stored at − 80 °C until sectioning. For sectioning, 10 μm consecutive sections were cut on an adhesive cryofilm type 3C (16UF) using a modified Kawamoto method 60 , to help preserve the morphology of the tissue during sectioning. The sections were left to dry for 10 min prior to use. The corneal sections and the primary limbal cells cultured on a chamber slide were fixed with 4% paraformaldehyde, permeabilized with 0.1% Triton X-100 in PBS for 10 min and blocked with 2% BSA solution in PBS anti-P63 (ΔNp63) (1:100 dilution; Abcam; ab735). The tissue sections and primary limbal cells were washed five times and incubated with secondary antibodies diluted in 2% BSA blocking solution, goat anti-rabbit A488 (1:300 dilution; Thermo Fisher Scientific) or goat anti-mouse A568 (1:300 dilution; Thermo Fisher Scientific), for 50 min at ambient temperature in the dark. Cell nuclei were stained with 0.5 μg/mL DAPI for 10 min. The samples were washed five times in PBS, mounted with coverslips with Fluoromount G mounting medium (Thermo Fisher Scientific) and examined on a Nikon Eclipse Ti-E inverted microscope equipped with a X-Light V2-TP spinning disk (Crest Optics).

Data availability
The scRNAseq dataset presented in this study has been uploaded to the Gene Expression Omnibus (GEO accession: GSE186433) and can be accessed using this link: https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE18 6433. All the datasets presented in this study are available in DataVerseNL and can be accessed using this link: https:// doi. org/ 10. 34894/ X7ZSDZ.