Single-Cell Transcriptomics Reveals Cellular Heterogeneity and Complex Cell–Cell Communication Networks in the Mouse Cornea

Abstract

Purpose: To generate a single-cell RNA-sequencing (scRNA-seq) map and construct cell–cell communication networks of mouse corneas.

Methods: C57BL/6 mouse corneas were dissociated to single cells and subjected to scRNA-seq. Cell populations were clustered and annotated for bioinformatic analysis using the R package “Seurat.” Differential expression patterns were validated and spatially mapped with whole-mount immunofluorescence staining. Global intercellular signaling networks were constructed using CellChat.

Results: Unbiased clustering of scRNA-seq transcriptomes of 14,732 cells from 40 corneas revealed 17 cell clusters of six major cell types: nine epithelial cell, three keratocyte, two corneal endothelial cell, and one each of immune cell, vascular endothelial cell, and fibroblast clusters. The nine epithelial cell subtypes included quiescent limbal stem cells, transit-amplifying cells, and differentiated cells from corneas and two minor conjunctival epithelial clusters. CellChat analysis provided an atlas of the complex intercellular signaling communications among all cell types.

Conclusions: We constructed a complete single-cell transcriptomic map and the complex signaling cross-talk among all cell types of the cornea, which can be used as a foundation atlas for further research on the cornea. This study also deepens the understanding of the cellular heterogeneity and heterotypic cell–cell interaction within corneas.

The transparent avascular cornea allows light into the eye and refracts light on to the retina. It also insulates internal ocular tissue from the external environment.¹ Structurally, it is composed of five layers: the outermost stratified epithelium, Bowman's layer, collagenous stroma, Descemet's membrane, and the innermost endothelium. The epithelium is an important barrier against external injury. The corneal stroma, composed of highly aligned collagen fibers and an extracellular matrix produced from embedded keratocytes, provides key mechanical support for the structural integrity of the cornea. The corneal endothelium consists of a monolayered hexagonal cell sheet residing on Descemet's membrane. It not only allows nutrient influx from the aqueous humor into the cornea but also actively pumps excess water from the corneal stroma back into the aqueous humor to prevent corneal edema.² 

Table.

 

The Detailed List of the Corresponding Antibodies and Desired Concentrations

 

Single-cell RNA sequencing (scRNA-seq) technology is a powerful tool that can reveal the heterogeneity of cells and their functions in complex biological systems of normal or diseased status.^3,4 Several human corneal single-cell transcriptomic maps have been reported.^5–9 The heterogeneity of cellular components of the human cornea has been revealed by scRNA-seq analysis of postmortem human cornea samples.⁹ These samples were often kept in culture/storage medium for different durations before sequencing, which might account for the variation in the reported transcriptomes. Mice are a popular model for corneal research owing to advances in genetic tools.¹⁰ A transcriptomic atlas of all corneal cell populations and a map of their cell–cell signaling networks could provide crucial information to analyze the physiological regulation and pathobiology of corneal diseases. In addition, such a data set can serve as a model map for future research on corneal physiology and pathology. Up-to-date scRNA transcriptomes of mouse corneal epithelium have been reported,^11,12 but the need remains for a complete transcriptomic atlas that includes and compares all cell populations within the mouse cornea. 

In this study, we constructed a complete scRNA-seq transcriptomic atlas of normal mouse corneas. We generated an scRNA-seq data set based on 14,735 corneal cells isolated in the physiological state. We identified 17 cell clusters, including all major cellular components within the cornea, delineated the cellular heterogeneity of the cornea, created individual spatial maps of cluster-specific marker expression using immunofluorescence (IF) analysis, and mapped the global intercellular signaling networks. Our data can be used as a reference for further study of corneal physiology and pathobiology. 

Material and Methods

Tissue Preparation

Eight-week-old C57BL/6 mice were purchased from the Taiwan National Laboratory Animal Center, Taiwan. All procedures with mice were approved by and performed in accordance with the requirements of the Institutional Animal Care and Use Committee of National Taiwan University. The study also adhered to the Association for Research in Vision and Ophthalmology Statement for the Use of Animals in Ophthalmic and Vision Research. 

Isolation of Single Cells From the Mouse Cornea

In the prior work on mouse single-cell study,¹¹ the EDTA-based digestion protocol yielded abundant corneal epithelial cells but a few keratocytes and very few, if not at all, corneal endothelial cells. To obtain a comprehensive transcriptomic map of cornea, we modified the cell isolation method to ensure that a sufficient number of all corneal cell types was obtained, especially endothelial cells and keratocytes. Detailed procedures can be found in the supplementary materials and methods section. After tissue dissociation, we used restricted gating cell sorting (FACSAria III; BD Bioscience, Franklin Lakes, NJ, USA) with FSC/SSC and PI-negative to exclude cell debris, doublets, aggregates and unqualified cells. 10x Genomics (Pleasanton, CA, USA) recommends cell sorting or the dead cell removal kit for eliminating cell debris to enhance cell viability. Cell sorting offers additional capability to remove cells between healthy and unhealthy. Thus we opted for cell sorting instead. A total of ∼1250 single-live keratocytes/corneal endothelial cells and ∼5250 single-live corneal epithelial cells were obtained from subsequent cell sorting. The cell number of corneal epithelium and corneal stromal cells/endothelium were then pooled at a ratio of 1:1 before loading on the chip. Although our cell viability decreased during long-time cell sorting, the cell viability before loading on the chip was 99%. The flowchart is summarized in Figure 1A. 

Figure 1.

 

ScRNA-seq analysis of adult mouse cornea. (A) Flowchart of the scRNA-seq experiment. (B) Transcriptomic profiles and unbiased clustering of two mouse cornea libraries shown by UMAPs for the 17 distinct clusters identified in this study. (C) Integrated UMAPs of the above two libraries because of identification of similar clusters. (D) Comparative heatmap of selected differentially expressed genes demonstrating the different transcriptomic profiles of the 17 cell clusters. (E) Feature plots demonstrating representative gene signatures of major cell clusters.

Single-Cell RNA Sequencing

Cells were captured and libraries were generated. This involved gel bead-in-emulsion generation, barcoding, post gel bead-in-emulsion–reverse transcription (RT) cleanup, cDNA amplification, library preparation, quality control, and sequencing with the use of Single Cell 3ʹV1 chemistry (10x Genomics). Library quality control metrics were as follows: 41.4 ng/µL with an average size of ∼495 bps. Libraries were sequenced to ∼40,000 reads per cell on the NovaSeq 6000 platform (Illumina, San Diego, CA, USA). All scRNA-seq data were deposited in the Sequence Read Archive database under the accession code PRJNA891516. All scRNA-seq data sets can be assessed and visualized in http://www.sjlin.tw. In addition, the detailed procedure for bioinformatics analysis of the scRNA-seq data, trajectory inference, cell-cycle discrimination, and cell–cell communication network analysis are provided in the Supplementary Material and Methods. 

Results

Single-Cell Transcriptomic Atlas of the Mouse Cornea Reveals Heterogeneity of Corneal Cells

We sequenced the transcriptomes of single corneal cells using the droplet-based scRNA-seq technique (Fig. 1A). We constructed two scRNA-seq data sets from 20 cornea samples each for analysis with the 10x Genomics platform. One data set yielded 8575 cells with 21,978 genes and the other 6157 cells with 25,991 genes. A relatively low cell number obtained resulted from additional cell sorting to achieve a 99% viability. Uniform manifold approximation and projection (UMAP) demonstrated similar transcriptomic profiles for the two data sets (Fig. 1B). Unsupervised clustering categorized the cells in both data sets into 17 transcriptionally distinctive clusters. 

We then integrated the two data sets for downstream analysis with high consistency (Fig. 1C). The distinguishable transcriptomic profiles of these 17 cell clusters are shown as representative heatmaps with differentially expressed genes (DEGs) (Fig. 1D). The 17 clusters were annotated using the gene expression pattern of known cell-type-specific marker genes as discussed below. The cell clusters corresponded to seven major distinctive subtypes according to their expression of known marker genes (Fig. 1E). Three cell clusters (1, 2, and 5) displaying high levels of keratocan (Kera), the specific marker of corneal stromal keratocytes (CSKs),¹³ were designated as CSK clusters. Two cell clusters (6 and 11) were identified as endothelial cells (CEnCs) because of their high expression of collagen 8A1 (Col8A1).¹⁴ One cluster (13) was identified as vascular endothelial cells because of its high expression of Pecam1 (CD31).¹⁵ One cluster was designated as immune cells because of F4/80 (Adgre1) expression.¹⁶ Another cluster exhibited a high level of Clec3b and Col14A1, suggestive of fibroblasts.¹⁷ 

The remaining large cell population expressing a moderate-to-high level of keratin 12 (Krt12) was assigned to the corneal epithelium.¹⁸ Residing close to Krt12⁺ cells were two small clusters (16 and 17) expressing the Krt13 gene. They were designated as conjunctival epithelial cells.¹⁹ The expression of Krt12, Kera, and Col8A1 represented conventional markers for three major corneal cell types (epithelium, keratocytes, and endothelium). We also performed trajectory analysis (Fig. 2A) and cell cycle analysis (Fig. 2B) to map the progression of differentiation and the cell cycle status, respectively. 

Figure 2.

 

Trajectory analysis of corneal cell subtypes. (A) Trajectory inference analysis of scRNA-seq transcriptomic profiles of 17 clusters. (B) Cell cycle scoring analysis of 17 cell clusters in the mouse cornea. (C) Distribution of root cells (left panel) and end points (right panel). (D) A representive pseudotime inference utilizing slingshot analysis with the starting point set as cluster 14 for corneal epithelial clusters.

Transcriptomic and Spatial Characterization of Epithelial Cell Subclusters

Nine clusters displaying moderate-to-high expression of Krt12 indicated their fate as the corneal epithelium (Fig. 3). To map the process of differentiation from limbal stem cells to terminally differentiated epithelial cells, we performed trajectory analysis of the nine epithelial cell clusters (Figs. 2C, 2D). We found expression patterns of epithelial subcluster marker genes that have been reported in previous mouse cornea scRNA-seq maps^11,12 and also found both distinctive and overlapping gene expression between epithelial subclusters (Fig. 3A). This finding agrees with human scRNA-seq data showing that epithelial cells may not only have exclusive markers but also share common gene expression as a part of the continuous epithelium differentiating from limbal-located stem cells to terminal differentiated cells.⁶ 

Figure 3.

 

ScRNA-Seq analysis of the corneolimbal epithelium. (A) Violin plots showing common and preferential expression of representative marker genes in epithelial clusters. (B–D) Immunofluorescence staining of the whole cornea to demonstrate differential expression of genes including (B) Ocln/Areg (C) MKi67/Krt14, (D) Knstrn/PBK, (E) Gpha2/PBK, (F) Krt12/PBK, and (G) Krt15. The direction of the yellow arrowhead indicates basal layer to superficial layer.

Two Krt12⁺ clusters (7 and 4) were identified as differentiated corneal epithelial cells (CEpiDCs) because of the high expression of tight junction-related genes (Cldn3 and Ocln) and suprabasal cell markers (Cldn4, Cdkn1a, and Dsg1a) and minimal expression of basal cell markers (Fig. 3A).¹¹ Both CEpiDC clusters withdrew from active cell cycling into a postmitotic state (Fig. 2B). CEpiDC II expressed higher superficial tight-junction genes, such as Ocln, and evolved from CEpiDC I in the epithelial trajectory (Fig. 2D). Therefore we assigned CEpiDC II to the terminally differentiated superficial epithelium and CEpiDC I to suprabasal cells. IF staining confirmed strong expression of Ocln within the superficial epithelium (Fig. 3B; Supplementary Video S1). 

To characterize the differentiation status of epithelial subclusters, we compared reported stem cell markers and cell cycle-related genes in the remaining five Krt12⁺ clusters (Fig. 3A). We first identified cluster 14 as limbal stem cells (LSCs) because of its high expression of the proposed limbal stem cell markers Gpha2, Ifitm3, and Krt15, low expression of proliferative markers including Mki67, Top2a, Ccnb1, Ccna2, Birc5, and Rrm2, high expression of basal cell markers of Krt14, Itgb1, and Itgb4.¹² Most cells in this cluster were arrested in a quiescent state whereas only a few cells were in the G2/M phase (Fig. 2B). Consistent trajectory results were obtained through both unsupervised and supervised analyses (Fig. 2), including an unsupervised scVelo analysis, as well as three additional analyses using slingshot to validate cluster 14 as the start point, the origin of cells that subsequently differentiated into three trajectories (Fig. 2D). Therefore cluster 14 was annotated as LSC, despite the low expression of the cornea-specific marker Krt12 and conjunctival marker Krt13, which was consistent with previous studies showing co-localization of differentiation and stem cell markers.^11,12 The IF staining confirmed the spatial localization of Gpha2 mostly at the basal layer of the outer limbus (Fig. 3F; Supplementary Video S2). 

Clusters 9 and 12, displaying high proliferative markers (Fig. 3A) with most cells in G2/M and S phases (Fig. 2B), were designated as highly proliferative transit-amplifying cells (TACs). The Krt12⁺ cluster 8 expressed a moderate level of proliferative markers, which was higher than their expression in CEpiDCs and LSCs but significantly lower than that in the two TAC clusters 9 and 12. Most cells in cluster 8 were in the S phase. According to trajectory analysis, cluster 8 stemmed from LSCs and developed into two trajectory branches. One branch led to cluster 9 (a highly proliferative TAC cluster), whereas the other led to the differentiated cells in clusters 3, 4, and 7, suggesting the inherent commitment of this branch to differentiation (Fig. 2D). Accordingly, we defined cluster 8 as an early TAC cluster (TAC I). One of the top DEGs of cluster 8 was the UNG gene, an important safeguard in the DNA repair system, which is highly expressed in human pluripotent stem cells.²⁰ IF staining showed a positive but somewhat stochastic UNG expression in corneas (Supplementary Fig. S1). 

The highly proliferative TAC cluster 12 stemming from LSCs expressed a moderate level of Krt15 and had a tendency to refuel the LSC cluster according to pseudotime data. We identified one of the top DEGs, Knstrn, in cluster 12 (Fig. 1D) and validated its spatial localization mainly, but not exclusively, in the inner limbal area (Fig. 3D; Supplementary Video S3). In addition, Krt12/Gpha2 IF co-staining revealed a Krt12⁻Gpha2⁻ cell population located in the inner limbus, between Krt12⁺Gpha2⁻ corneal epithelial cells and Krt12⁻Gpha2⁺ LSCs, similar to the spatial expression of Knstrn+ cells (Fig. 1F). We identified cluster 12 as a highly proliferative TAC cluster located mainly in the inner limbus and directly originating from LSCs (TAC II). 

PBK has been defined as a mature TAC marker in the cornea, but not the limbus, in humans.¹¹ Originating from an early TAC cluster (cluster 8), the highly proliferative TAC population (cluster 9) exclusively expressed a high level of PBK; it was designated as a mature TAC population in corneas (TAC III). IF staining showed PBK expression in the inner limbus/ peripheral cornea area, but not in the central cornea. Notably, PBK was distributed more centripetally than Knstrn in cluster 12 (Figs. 3D, 3E; Supplementary Videos S3, S4). 

The last Krt12⁺ cluster (3), stemming from cluster 8 in the trajectory, expressed low levels of differentiation (Ocln, Cldn3, and Areg), stem cell (Gpha2 and Krt15), and proliferation (Mki67, Top2a, Ccnb1, Ccna2, Birc5 and Rrm2) markers and high levels of basal cell markers (Gja1 and Krt14). This cluster was designated as corneal basal epithelial cells (CEpiBCs) (Fig. 3C; Supplementary Video S5). Two small clusters (16 and 17) embedded within Krt12⁺ corneal epithelial clusters were designated as conjunctival epithelial cells because of the expression of known conjunctival markers Krt13, representing some inevitable contamination of conjunctival epithelium in our specimens. 

Overall, by integrating whole-mount IF staining and transcriptomic profiles, we identified and spatially mapped subtypes of epithelial lineages from limbal stem cells to terminally differentiated cells, including one LSC cluster, one low-proliferative early TAC cluster, two highly proliferative TAC clusters, and three differentiated cell clusters. We also found that LSCs produced dissimilar differentiation branches (Fig. 2D). 

Transcriptomic and Spatial Characterization of Keratocytes

We identified three clusters (1, 2, and 5) expressing high keratin (Kera) as CSKs (Fig. 1C). The Kera gene is critical for maintaining corneal structure.¹³ Trajectory analysis showed the differentiation progression from CSK I to CSK III (Fig. 2A). Cell cycle analysis showed all three keratocyte clusters at G1 phase or S phase (>90%) (Fig. 2B). Next, we compared the representative DEG expression of these three subtypes and found that Egfl6 was enriched in CSK I (cluster 1) whereas Col14A1 was enriched in CSK II and CSK III (clusters 2 and 5) (Fig. 4A). Whole-mount IF staining showed the presence of Col14A1 across the whole stroma, albeit higher in the posterior stroma, whereas the ECM protein Egfl6 was preferentially expressed in the anterior stroma (Fig. 4B; Supplementary Videos S6, S7). In addition, heatmap analysis showed fibroblasts with high Col14A1 (Fig. 1D) and IF staining (Supplementary Video S8). 

Figure 4.

 

ScRNA-Seq analysis of keratocytes and the corneal endothelium. (A) Violin plots showing common and preferential expression of representative marker genes in keratocyte clusters. (B) Immunofluorescence staining of the whole cornea to demonstrate differential expression of genes, including Col14A1 and Egfl6, indicating different distributions within the stroma. (C) Violin plots showing common and preferential expression of representative marker genes in endothelial clusters. (D) Immunofluorescence staining of the whole cornea revealing the heterogenous expression of CD109 in endothelial cells.

Transcriptomic and Spatial Characterization of Corneal Endothelial Cells

ScRNA-seq information for the mouse corneal endothelium is lacking because the previously reported scRNA-seq transcriptomic profiles of the corneal endothelium are from human samples.^6,21 Unsupervised clustering identified two cell clusters (6 and11) as CEnCs with high expression of Col8A1, rather than Col8A2 as in human corneas (Fig. 1D).¹⁴ Col8A1 as a CEnC marker was validated by its exclusive localization in the corneal endothelium by IF analysis (Fig. 1F). 

The two CEnC clusters were partially dispersive in their transcriptional profiles (Fig. 4C). Both clusters were mostly arrested in the G0/G1 phase whereas a small proportion was in the S phase (Fig. 2B). This was consistent with the limited proliferative capacity of the endothelium in vivo.²² The CEnC I cluster showed higher expression of mitotic inhibition–related genes (Cdkn1a, some TGF-βs, and CD109), whereas CEnC II exhibited higher expression of endothelial–mesenchymal-transition (EnMT)–related genes (Fndc1, Acta2, and Myoc). The transcriptome differences suggest functional heterogeneity in the corneal endothelium. IF staining of candidate DEGs, including Penk, Chad, Prss12, Smoc1, and Postn, was unable to spatially differentiate between the two clusters. We then examined the expression of EnMT-related CD109 and found that CD109^high cells were scattered in the corneal endothelium in a somewhat stochastic pattern (Fig. 4D; Supplementary Video S9). This observation suggested the presence of cellular heterogeneity in the corneal endothelium. 

Transcriptomic and Spatial Characterization of Vascular Endothelial Cells, Immune Cells, and Fibroblasts

Within the 17 cell clusters, we also identified three supportive cell clusters. Vascular endothelial cells were identified by the conventional marker Pecam1 (CD31) (Fig. 5A), with high expression of other characteristic markers, including Cldn5 and Col4A2 (Supplementary Video S10). Imaging of the whole cornea showed distribution of vascular endothelial cells at the periphery of corneolimbal samples (Fig. 5B; Supplementary Video S11), representing limbal vasculature. We identified immune cells by Adgre1 (F4/80) gene expression, a monocyte/macrophage marker (Fig. 5C), which mainly accounts for innate immunity of the cornea.²³ This immune cell cluster presented a mixed immune cell phenotype of Langerhans cell markers (Ccl3, Msr1, and Ptprc)²⁴ and a macrophage-specific marker (Fcgr1).²⁵ The distribution of F4/80⁺ immune cells decreased centripetally as shown in IF images of the whole cornea (Fig. 5D; Supplementary Video S12). Langerhans cells, with their characteristic dendritic morphology, reside in the basal epithelial layer of the peripheral corneal and limbal area. The corneal stroma and the basal layer of the corneal epithelium were enriched with F4/80⁺ macrophages with a pleomorphic morphology (Fig. 5D). We also characterized fibroblasts with Col14A1 and Clec3b expression based on mouse scRNA-seq data²⁶ (Fig. 5E). Whole-cornea IF imaging showed Clec3b⁺ cells in the limbal stroma and corneal endothelium (Fig. 5F; Supplementary Video S8). 

Figure 5.

 

ScRNA-Seq analysis of vascular endothelial cells, immune cells, and fibroblasts. Violin plots of representative differentially expressed genes for three minor cell types in the cornea: (A) Pecam1 (CD31) for vascular endothelial cells; (C) Adrge1 (F4/80) for immune cells; and (E) Clec3b for fibroblasts. (B) Whole-cornea immunofluorescence staining for CD31. White and yellow arrowheads indicate the corneal epithelium and corneal endothelial cells, respectively. (D) Whole-cornea immunofluorescence staining for F4/80. (F) Whole-cornea immunofluorescence staining for Clec3b. White and yellow arrowheads indicate fibroblasts and corneal endothelial cells, respectively.

CellChat Inference of Cell–Cell Communication Between Multiple Cell Populations in the Mouse Cornea

We used CellChat²⁷ to study the global intercellular communication networks within mouse corneas (Fig. 6). We first constituted and quantified the global signaling cross-talk atlases, which provided the number of cell–cell interactions and their strengths in circle plots (Figs. 6A, 6B). The thicker the line was, the higher the number of interactions/strength. There were complex cell–cell interactions within the cornea (Figs. 6A, 6B), with autocrine signaling loops present in all cell types. In the global view, we found that CSK I, CSK II, and CEpiBCs accounted for most intercellular communications (Fig. 6A). The outgoing signaling between 39 most significant pathways in CellChat was studied using pattern recognition analysis (Fig. 6C). We chose one major signaling pathway with high communication strength for the epithelium, keratocytes, and endothelium, respectively. 

Figure 6.

 

CellChat-inferred cell–cell communication networks reveal the functional heterogeneity of corneal cell populations. Global communications are presented by circle plots showing the number (A) and weight/strength (B) of significant ligand–receptor pairs in 17 cell clusters. The edge width is proportional to the indicated number/weight of individual ligand–receptor pair signaling. The loop in the individual plot indicates the autocrine pathway of the cell cluster. (C) The dot plot shows the comparison of outgoing signaling patterns in top 39 signaling pathways of secreting cells in the 17 cell clusters. The dot size is proportional to the contribution score obtained from pattern recognition analysis. A higher contribution score indicates that the signaling pathway is more enriched in the corresponding cell cluster. Three representative signaling pathways, EGF (D–F), ncWNT (G–I), and TGF-β (J–L) pathways, were further analyzed. The inferred networks of communication between all cell types are displayed using hierarchical plots (D, G, J). The relative contribution of each ligand–receptor pair to the global signaling networks (E, H, K), and violin plots showing the expression distribution of signaling genes in each signaling network (F, I, L) are also shown.

Network centrality analysis of the inferred epidermal growth factor (EGF) signaling network showed that epithelial cells were the most prominent EGF ligand sources and were themselves the primary EGF signaling targets (Fig. 6D). These findings are consistent with the known critical role of EGF in corneal epithelial homeostasis.²⁸ The relative contribution of ligand–receptor pair analysis showed that the HbEGf–EGFR pair contributed most to the EGF signaling pathway (Fig. 6E), which is consistent with a previous report.²⁹ ScRNA-seq violin plots showed that EGF signaling gene expression dominated in epithelial clusters (Fig. 6F). Notably, LSCs were not a source of EGF ligands but only the target of EGF signaling. This finding confirms that EGFR inhibition affects LSCs proliferation and further corneal stratification during corneal epithelial wound healing and also may play a role in maintaining the normal thickness of the epithelium.³⁰ 

CellChat-inferred ncWNT signaling network revealed that keratocytes and some endothelial cells (CEnC I) were the main ligand sources (Fig. 6C), which acted both autocrinally and paracrinally (Fig. 6G). WNT5a–Fzd2 was the major ligand–receptor pair contributing to the ncWNT communication networks (Fig. 6H). ScRNA-seq violin plots showed that ncWNT-pathway-related gene expression dominated in keratocyte clusters and that expression was minor in the endothelium (Fig. 6I). As in the skin,²⁷ the FGF signaling network was similar to that of ncWNT (Supplementary Fig. S2A) with a predominant FGF18–FGFr1 ligand–receptor pair (Supplementary Fig. S2B). This result agrees with previous work showing that FGF-1 and -2 promote the fibroblast phenotype and reverse the myofibroblast phenotype induced by TGF-β in keratocytes.³¹ 

TGF-βs play a key role in corneal development, healing of corneal scaring,³² and maintenance of the corneal endothelium in a quiescent state.³³ TGF-β proteins are expressed in the corneal epithelium and are also present in the aqueous humor.^32,34,35 Our CellChat data showed that corneal endothelial cells were the main TGF-β source acting autocrinally and paracrinally (Fig. 6J). CellChat-inferred TGF-β networks showed a discrepancy in the roles of the epithelial lineage participating in the TGF-β pathway (Fig. 6J). Terminally differentiated superficial epithelial cells did not participate in the TGF-β pathway. The LSCs acted as receivers only for TGF-β signals. Epithelial subtypes, except for CEpiDC II and LSCs, were TGF-β ligand sources, wheras only TAC III and LSCs acted as receivers. Keratocytes were both sources and targets in the TGF-β pathway and communicated with both the epithelium and the endothelium, which confirmed its role in regulating fibrosis in wound healing. CEnC clusters were both autocrinal and paracrinal targets and sources (Fig. 6J). CEnC II contributed more as a TGF-β source compared to the CEnC I cluster. In addition to autocrine, CEnCs received TGF-β signaling from most corneal cell types except for LSCs and CEpiDC II. The major ligand–receptor in TGF-β signaling networks was the activin A receptor like type 1 (Acvrl1)–TGF-βr1 pair (Fig. 6K). Violin plots showed that, in addition to minor clusters of fibroblasts and vascular endothelial cells, CEnC clusters had higher expression of all TGF-β and TGF-βr isoforms (Fig. 6L). The TGF-β1r gene was expressed in the epithelium, keratocytes, and endothelium; this was consistent with the known role of TGF-β signaling in corneas.^32,33 

CellChat revealed four novel significant pathways in the healthy mouse cornea, including Granulin precursor (Grn), Pros1(Pros), chemerin, and visfatin. Grn binds to sortilin (Sort1) to regulate lysosomal function and microglial responses in the CNS.³⁶ Our results showed Sort1 expression in 3 keratocytes clusters and CEpiDCs II, whereas progranulin (Pgrn) was expressed in all cell clusters (Supplementary Figs. S3A–C). Pros1 and its receptor AxL receptor tyrosine kinase (Axl) pathway regulate intrinsic mesenchymal signaling and the extrinsic immune microenvironment, contributing to the growth of aggressive tumors.³⁷ Our results revealed the autocrine manners in three keratocyte clusters, two corneal endothelial cell clusters, and one fibroblast cluster. The Pros1 is expressed in the vascular endothelial cell cluster, and the Axl is expressed in immune cells (Supplementary Figs. S3D–F). 

Discussion

Unlike the averaged global gene expression determined by bulk RNA sequencing, scRNA-seq allows precise and simultaneous mapping of gene expression patterns in individual cells.³⁸ We provide a comprehensive scRNA-seq atlas of the healthy mouse cornea, with 17 distinct cell clusters identified and spatially mapped, and unveil the complex cell–cell communication networks using CellChat. We successfully designated the corneal endothelium and the small population of LSCs. No previous mouse scRNA-seq analysis has reported comprehensive endothelial transcriptomic profiles.^5–8,11,12 

Several scRNA-seq analyses of human and mouse corneal epithelial lineages, including our study, have disclosed an undefined differentiation program, with a certain degree of both interspecific and intraspecific discrepancies.^{5–8,11,12,39} Differences between species, diversity in scRNA-seq methodologies, and the lack of confirmed corneal epithelial staged markers might have led to the variable results. We identified seven clusters of corneolimbal epithelial lineages in normal mouse cornea based on the expression of known maturation, proliferation, and reported but inconclusive, stem cell markers, and validated their differentiation status using trajectory analysis (Fig. 2). In our study, only one cluster met LSC criteria, including low corneal and conjunctival maturation markers, low proliferative markers, basal localization, and positive proposed putative LSC marker gene expression (Gpha2⁺/Ifitm3⁺/CD63⁺/Krt15^high/Krt14⁺/Krt12^low/Muc20⁻) (Fig. 3A). The LSC cluster also expressed Id1/3 and Sfrp1¹² genes, compatible with LSC characteristics. Not all cells of the LSC cluster were at the quiescent G1 phase (Fig. 2A), compatible with Altshuler's data and the equipotent stem cell model hypothesis that LSCs are not only found at G0/G1 phase but also at S phase.¹² Although we identified a LSC cluster that was compatible with the “quiescent LSC” population in the bi-compartmented LSC theory,¹² we did not identify an “active LSC” cluster exclusively expressing the proposed active LSC markers.¹² Instead, we identified a TAC cluster (TAC I) stemming from the LSC cluster with Mki67^high40/Birc5^high41/Krt15^mod/Atf3^high/Mt1^high expression (Fig. 3A) located in the inner limbal area. Therefore we inferred that the TAC I cluster might overlap with the “aLSC” population in Altshuler's hypothesis,¹² despite the discrepancy in DEG expression. In trajectory analysis, the TAC I cluster was inferred not only to stem from LSCs but also to fuel the LSC cluster (Figs. 2A, 2D). Our inference was not completely consistent with the bi-compartmented LSC model.¹² Therefore further studies are needed to clarify this inference. In Altshuler's study, two “mitotic cell” clusters composed of mixed cell types, but not TACs, were annotated.¹² In our work, we identified three TAC clusters based on their moderate-to-high proliferative marker gene expression (Fig. 3A) and inferred their differentiation stages by trajectory analysis. We also identified differentiated corneal basal, suprabasal, and terminally differentiated superficial epithelial cells (Fig. 3). We did not perform subclustering of superficial/suprabasal/basal limbal epithelium as in previous studies with human data,^5,7,8 because of the small number of cells, depth of sequencing, and the sequencing approach. 

We identified three clusters of keratocytes (Fig. 1C) and revealed differences in DEG expression between the clusters: Egfl6 expression was enriched in CSK I with preferential anterior stromal distribution (Fig. 4B), whereas Col14A1 expression was enriched in CSK II and CSK III with the protein distributed across the full thickness of the stroma, especially the posterior stroma (Fig. 4A and B). Egfl6 is known to be expressed in the early stage of development and is related to tumorigenesis and angiogenesis.⁴² It implies the contribution of CSK I to the development of the stroma and a possible correlation with angiogenesis following injury.⁴² Collagen XIV is crucial for regulating fibrillogenesis and tissue mechanics in the corneal stroma during development and following injury and is also responsible for the maturation of the endothelium.^43,44 The preferential expression of the Col14A1 gene in CSK II and CSK III reflects their roles in regulating stromal architecture and biomechanical strength. We did not annotate a subpopulation of corneal stromal stem cells (CSSCs) with Kera^low/Fhl1^high/Thy1^high/Six2⁺/Scf⁺ expression as in a previous study of mice,¹¹ nor with MMP3^high expression as in a study of human data.⁵ In contrast to Fhl1^low/Thy1^low keratocytes, Clec3b⁺ fibroblasts showed significantly upregulated expression of Fhl1, Thy1, and Fbln1 (human limbal fibroblast marker).⁵ Therefore we designated the Clec3b^high/Fhl1^high/Thy1^high/Fbln1^high/Kera^low/ Six2⁻/Scf⁻ cluster as limbal fibroblasts, not as CSSCs. In our study, the expression of potential CSSC markers was similar between keratocyte clusters in UMAP. However, we identified the same novel keratocyte DEG of NNMT⁶ as a common keratocyte marker for human and mice corneas. 

We provide the first mouse corneal endothelial transcriptomic profiles. Most reported human endothelial marker genes^6,9 did not show significant expression in the mouse endothelium, whereas TSPAN6, FZD2, and FGFR1 expressed moderately, but not exclusively, in the mouse cornea. Unlike the human endothelium⁹ where many respiratory markers showed significantly higher expression than in other cell types, only ATP1a1 and SLC4A4 expressed significantly, but not exclusively, in the mouse endothelium. Similar to the human endothelium,⁶ both endothelial clusters expressed the PITX2 gene, whereas CEnC II showed higher expression than CEnC I. The expression of EnMT-related genes varied markedly in the human endothelium,^5,6,9 possibly because of sample preparation.⁶ In our study, the expression of EnMT-related genes (ACTA2 and TAGLN) was higher in CEnC II than in CEnC I at low-to-moderate levels. We suggest that our endothelium data reflect the normal condition without the influence of sample preparation⁵ and that the CEnC II cluster contributed more to the EnMT process. CEnC I expressed higher mitotic inhibition-related genes (Cdkn1a, some TGF-βs, and CD109), reflecting its role in the maintenance of endothelial homeostasis. Both mouse keratocyte DEGs and endothelial subpopulations were different from human data.⁶ Nevertheless, we demonstrated the first mouse keratocytes and endothelial transcriptomic profiles with the presence of functional heterogeneity. 

Another major outcome of this study is that we demonstrated the complexity of signaling between corneal cells. We provided the first CellChat analytic report on mice cornea, revealing quantitative inferences of intercellular signaling and predicts inputs and outputs.²⁷ We displayed the top 2% ranking of strong cell-cell interactions in cornea, including well-known pathways and some novel pathways that have not been addressed before, such as Pros, Grn, chemerin and visfatin. (Fig. 6C and Supplementary Fig. S3). First, we studied EGF signaling networks. All three major corneal cell types may express EGF and its receptor,^30,45,46 with EGF receptors also present in limbal basal cells.⁴⁷ In the cornea, EGF signaling plays crucial roles in maintaining epithelial homeostasis,²⁸ myofibroblastic differentiation and migration,⁴⁸ and endothelial wound healing, and as a mitogen for the corneal endothelium.⁴⁹ According to the CellChat analysis, the epithelium was the major source and target of EGF ligands, mostly autocrinally. Of note, LSCs receive EGF signals only, confirming their role as a repository of cells for regeneration of the epithelium. EGF signaling was previously thought to act partially and autocrinally in the three major corneal cell types because of the coexistence of EGF ligands and receptors.^48,50 However, network centrality analysis decrypted that in fact, subpopulations of keratocytes and endothelium act differently in EGF signaling networks (Fig. 6D). 

Compared to the significant role of the canonical WNT pathway in regulating LSC proliferation,⁵¹ the ncWNT pathway plays roles in migration and development of the endothelium and the pathophysiology of keratoconus, although the detailed mechanism remains unclear.^52–54 Our CellChat data inferred that keratocytes were the main ncWNT sources both autocrinally and paracrinally (Fig. 6). Similar to the communication between dermal fibroblasts,²⁷ we found cell–cell communication networks of FGF and the ncWNT pathway (Supplementary Fig. S2A). These results suggest different roles for canonical WNT and ncWNT signals in corneal homeostasis. 

The TGF-β family plays a complex role in all three major cell types in the cornea for maintaining corneal integrity and in the wound healing process.⁵⁵ TGF-β has been shown to be restricted in the epithelium of the normal healthy cornea. TGF-β1 is responsible for myofibroblastic differentiation during stromal wound healing.⁵⁵ TGF-β2 in the aqueous humor is vital for inhibiting corneal endothelial proliferation by suppressing entry into the S phase during homeostasis.^56,57 We confirmed that all three TGF-β isoforms participated in corneal physiology (Fig. 6K). CellChat inferred that functionally heterogeneous epithelial subpopulations played different roles in the TGF-β pathway (Fig. 6J), which has not been addressed before. Keratocytes were also main targets, reflecting their role in regulating fibrosis. The endothelium, but not the epithelium, acts as the major cells secreting TGF-β ligands in the normal cornea (Fig. 6C). The endothelial clusters are functionally diverse (Fig. 6J). CEnC II contributed more as TGF-β sources than the CEnC I cluster, consistent with the scRNA-seq finding that CEnC II expresses more EnMT-related genes. Although the exact nature of the newly discovered Pros, Grn and visfatin signaling pathways in cornea remains unclear, it has been observed that these three signaling pathways were associated with retina, including that Pros1 mutation correlates with inherited retinal degenerations,⁵⁸ Grn deficiency results in photoreceptor degeneration in mouse retina⁵⁹ and visfatin expression in diabetic rat retina.⁶⁰ Although chemerin has been shown to promote vessel growth in the mouse corneal angiogenesis model,⁶¹ its specific role in the cornea remains unexplored. 

In brief, the CellChat analysis provides a comprehensive map of the global cell–cell communication networks within the cornea. A revisit of the known cell-cell crosstalks may help to decrypt the signaling pathways in cornea and also functional heterogeneity in subpopulations of corneal cells. 

Conclusions

We provide a comprehensive scRNA-seq transcriptomic atlas of all cell types in the mouse cornea, and also reveal the complex intercellular signaling networks in the cornea using CellChat analysis. These bioinformatics data provide a fundamental blueprint for further research on the cornea and improvement of cell-based regenerative medicine for curing corneal blindness. 

Acknowledgments

The authors thank the Flow Cytometric Analyzing and Sorting Core and the Sequencing Core, Department of Medical Research, National Taiwan University Hospital, for technical support. 

Supported by Taiwan Bio-Development Foundation (TBF; to S.J. Lin), Chang Gung Memorial Hospital (CMRPG3G1623 to H.Y. Tan), Taiwan National Science and Technology Council (109-2314-B-182A-050 and 110-2314-B-182A-116-MY2 to H.Y. Tan; 111-2327-B-002-015 to S.J. Lin) and National Taiwan University Hospital (111IF0006; to S.J. Lin). Y.F. Wu was supported by a postdoctoral fellowship from Taiwan National Science and Technology Council (111-2811-B-002-156). This work was supported by the Brain Research Center under the Higher Education Sprout Project, co-funded by Taiwan Ministry of Education and the Ministry of Science and Taiwan National Science and Technology Council. The authors also thank the members of the S.J. Lin laboratory for their discussion. 

Disclosure: Y.-F. Wu, None; N.-W. Chang, None; L.-A. Chu, None; H.-Y. Liu, None; Y.-X. Zhou, None; Y.-L. Pai, None; Y.-S. Yu, None; C.-H. Kuan, None; Y.-C. Wu, None; S.-J. Lin, None; H.-Y. Tan, None