The Spatial Transcriptomic Atlas of Human Limbus and Vital Niche Microenvironment Regulating the Fate of Limbal Epithelial Stem Cells

Abstract

Purpose: This study aimed to generate the spatial atlas of the human limbus using spatial transcriptomic technology and reveal the deep interaction among the niche microenvironment.

Methods: The spatial transcriptomic atlas of human limbus was performed using 10× Genomics Space Ranger software platform. Single-cell RNA sequencing data of human limbal epithelial stem cells (LESCs) were downloaded for integrating analysis.

Results: We profiled more than 400 spots within each sample and spatially located major cell types within the limbus area. LESCs were localized mainly in the basement membrane, and limbal niche cells were situated predominantly within the stromal area. Next, the limbus was divided into four regions based on histological structure, and the differential expressed genes among the four regions were analyzed. Notably, GPHB5 was highly expressed in the epithelium of the middle region and co-staining with deltaNp63 suggested it might be a novel potential biomarker of LESCs. Subsequently, limbal mesenchymal stem cells were found to exhibit the greatest amounts of ligands associated with LESCs. The widespread activity of COL6A2/CD44 signaling among limbal mesenchymal stem cells, melanocytes, immune cells, and LESCs indicate its essential role in mediating bidirectional communication via the collagen pathway.

Conclusions: This research mapped the spatial positioning of key cells within the limbal niche and detailed interactions between major cell types. These findings provide a foundation for further LESC research and enhance our understanding of corneal biology.

The cornea is an avascular, transparent tissue that plays an important role in maintaining vision and protecting the eye from damage.¹ The homeostasis and renewal of the corneal epithelium are maintained by limbal epithelial stem cells (LESCs), which are uniquely localized to the Vogt palisade in the limbus and possess the key characteristics of epithelial stem cells, namely, undifferentiated, static, self-renewal, high proliferative potential, and tissue regeneration.² 

Like stem cells in other tissues, the functional activity of LESCs is regulated tightly by the surrounding microenvironment, also known as the niche.³ The niche is not only a protective environment, but also a necessary factor to maintain the characteristics of LESCs and initiate differentiation pathways in response to external signals, which determines cell survival, proliferation and even reverse differentiation.⁴ The main components of limbal niche include the extracellular matrix (ECM), chemical molecules, and limbal essential cells, such as mesenchymal stem cells (MSCs), melanocytes, and immune cells.⁵ MSCs interact with LESCs via various chemical and signal transduction pathways, as well as intercellular contacts and paracrine secretion.⁶ The main function of melanocytes is to transfer melanosomes containing melanin to LESCs to protect them from ultraviolet radiation and free radical damage.⁷ Immune cells, such as T cells, play important roles in maintaining LESC quiescence and corneal wound healing.⁸ Their interactions maintain metabolic homeostasis of the limbus and provide dynamic support for the proliferation, migration, stemness maintenance, and differentiation of LESCs.² 

In recent years, single-cell sequencing has dissected cell composition and cell–cell interaction networks at the single-cell level, opening a new window to reveal the heterogeneity and molecular regulatory properties of LESCs in the limbal basal epithelium, as well as new cell types in their progeny.^9,10 The technology of spatial transcriptomics retains tissue structure information to conduct comprehensive analysis of cells, such as in spatial organization of the mouse retina.¹¹ Therefore, to deeply understand the complicated limbal niche microenvironment and the underlying mechanism in regulating the LESCs fate, we applied the spatial transcriptomic technology to generate the spatial atlas of the human limbus. By integrating with reference single-cell RNA sequencing (scRNA-seq) data, we identified and analyzed major cell types in the limbus, as well as their spatial distribution and unique ligand–receptor connections. Collectively, our research provides a highly integrated picture of human limbus and transforms the understanding of LESCs at spatial level, which will help us to gain a deeper understanding of the niche microenvironment. 

Materials and Methods

Tissue Preparation

The collection of human tissues was approved by the Eye Bank and Medical Ethics Committee of the Ninth People's Hospital, Shanghai Jiao Tong University School of Medicine (SH9H-2021-T171-1), and was performed in strict accordance with the approved guidelines. After removal of the central buckle for corneal transplantation, two donor limbal rings (51 years old; cause of death, sudden death; time of absence, within 24 hours; and 61 years old; cause of death, car accident; time of absence, within 24 hours) were collected, and the conjunctiva and outer sclera at the limbal margin were removed under a microscope. 

10× Visium Spatial RNA-seq Data Preprocessing

Freshly collected human limbal tissues (Supplementary Fig. S1) were divided into appropriate size and embedded in Optimal Cutting Temperature and quickly frozen on dry ice. Tissue sections were subjected to methanol fixation, HE staining, imaging and destaining following the 10× Genomics recommended experimental procedure (CG000614). According to the 10× Genomics experimental flow (CG000495), probe hybridization, probe release, and transferred to the 10× Visium CytAssist slide, and the library construction was performed using the Visium CytAssist Spatial Gene Expression for FFPE kit. The DNA libraries were subjected to high-throughput sequencing using the PE-150 mode. 

By using 10× Genomics' Space Ranger software (version 2.0.1), the FASTQ files were processed and aligned to the GRCh38 human reference genome. Then, unique molecular identification counts were summarized for each barcode. The tissue overlay points were detected based on the images to distinguish them from the background. Subsequently, the filtered unique molecular identification count matrix was analyzed using the Seurat (version 4.1.0) R package.¹² 

The data were normalized using Sctransform¹³ to identify the top 3000 highly variable genes. Principal component analysis was used for reducing the dimension of the top-level variable gene–barcode matrices after logarithmic transformation. Graph-based clustering was conducted to cluster cells according to their gene expression profile with the FindClusters function. The two-dimensional uniform manifold approximation and projection (UMAP) algorithm with RunUMAP function was used to visualize the cells. The FindAllMarkers function (test. use = bimod) was performed to identify the marker genes in each cluster. 

Differentially expressed genes (DEGs) were screened using the FindMarkers function (test. Use=presto). An adjusted P value of less than 0.05 and a |log2foldchange| of greater than 0.58 were used as the thresholds of significantly differential expression. Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of DEGs were respectively conducted using R (version 4.0.3) based on the hypergeometric distribution. 

The Division of Four Regions

The division of the four regions is based on the results of hematoxylin and eosin staining. The region within the limbal epithelium is defined as the inner region. The region outside the limbal epithelium is defined as the outer region. In addition, based on the boundary between limbal epithelium and mechanism layer, the upper part is defined as epithelial layer of middle region, and the lower part is defined as stroma layer of the middle region. 

Robust Cell Type Decomposition (RCTD) Analysis

To infer the cell-type composition of each spot, we applied RCTD (version 1.1.0).¹⁴ When running RCTD, the creat.RCTD function used the default parameters, but each cell type had at least one cell and each spot had at least one unique molecular identifier, and the doublet mode parameter was set to FALSE at run time. Thus, the composition of cell types at each spot site was inferred. 

Cell-Cycle Analysis

The cell-cycle phase of individual cells was predicted using Seurat (version 3.1.2) with the CellCycleScoring function.¹⁵ Briefly, genetic markers for G2/M and S phase were used to assign cell scores, and cells expressing neither of the G2/M nor S phase markers were categorized as being in G1 phase. 

Monocle2 Pseudotime Analysis

Using Monocle2 (version 2.9.0) package with developmental pseudotime, specific steps are as follows.¹⁶ First, the importCDS function of Monocle2 package was used to convert the Seurat object into CellDataSet object, and the differentialGeneTest function of the Monocle2 package was used to select ordering genes (qval < 0.01), which were likely to be informative in the ordering of cells along the pseudotime trajectory. Then the reduceDimension function was performed for clustering, followed by trajectory inference with the orderCells function using default parameters. 

Identification of Transcription Factors (TFs) Using Single-cell Regulatory Network Inference and Clustering (SCENIC)

SCENIC analysis was performed using RcisTarget's motifs database and GRNboost (SCENIC version 1.2.4, RcisTarget version 1.10.0 and AUCell version 1.12.0) were run with default parameters.¹⁷ Combined with the base sequence in gene, TFs were identified in the list. The activity of each group of regulons in each cell was scored by the AUCell package (version 1.12.0). To evaluate the cell type specificity of each predicted regulon, the regulon specificity score, which is based on the Jensen–Shannon divergence, a measure of the similarity between two probability distributions was calculated. To be specific, the Jensen–Shannon divergence between each vector of binary regulon activity overlaps with the assignment of cells to a specific cell type was calculated. The connection specificity index for all regulons was calculated with the scFunctions package (https://github.com/FloWuenne/scFunctions/). 

Cell–Cell Communication Analysis

We used the CellPhoneDB (version 4.1.0)¹⁸ to identify biologically relevant ligand-receptor interactions from single-cell transcriptomics (scRNA-seq) data. A ligand or a receptor was defined as “expressed” in at least 10% of cells of a particular cell type and with a P value of less than 0.05 were calculated and analyzed. R packages Igraph and Circlize were used to display the cell–cell communication networks. 

The CellChat (version 1.1.3) R package¹⁹ was also used to analyze the cell communication. First, the CellChat object with the create CellChat function was created by importing the normalized expression matrix. The data were then preprocessed with identifyOverExpressedGenes, identifyOverExpressedInteractions, and projectData functions using the default parameters. The computeCommunProb, filterCommunication (min. cells = 10), and computeCommunProbPathway functions were then conducted to determine any potential ligand–receptor interactions. At last, the cell communication network was aggregated using the aggregateNet function. 

Immunofluorescence Staining

Limbal tissues were fixed with 4% paraformaldehyde, embedded in paraffin, and serially sectioned. The sections were incubated with 5% goat serum for 1 hour at room temperature and then treated with ALCAM (Proteintech, Rosemont, IL, USA; 21972-1-AP), PERP (Signalway Antibody, Greenbelt, MD, USA; NO. 49745), GPHB5 (Signalway Antibody, NO. 48970), ΔNp63 (Cell Signaling Technology, Danvers, MA, USA; 67825), CK14 (Abcam, Cambridge, UK; ab119695), CK3 (Santa Cruz Biotechnology, Dallas, TX, USA; sc-80000), CK12 (Proteintech, 24789-1-AP), MELANA (Proteintech, 60348-1-Ig), and TYRP1 (HUABIO, Taipei, Taiwan; ET7106-66) overnight at 4°C. Subsequently, the Alexa Fluor–conjugated secondary antibodies were stained for 1 hour at 37°C. Nuclei were stained with DAPI for 20 minutes and pictures were taken by fluorescent microscope. 

Transfection

One day before transfection, an appropriate number of LESCs were inoculated into the cell culture plate, and medium without antibiotics was added to each well to achieve a cell density of approximately 50% at the time of transfection. LESCs were transfected using Lipofectamine 3000 with 50 nM silencing RNA targeting GPHB5 according to the manufacturers’ instructions (Helix Biotech, Inc., Knoxville, TN, USA). LESCs were collected 72 hours after transfection. 

Statistical Analysis

In addition to bioinformatics analysis of the scRNA-seq data, all biological experiments were compared between the two groups using the Student t test. Three or more groups were evaluated using ANOVA tests, and the group mean pairs were compared with appropriate posterior tests. 

Data and Code Availability

All raw data are available in the Gene Expression Omnibus database with accession number GSE262250. 

Results

Spatial Transcriptomics Combined With scRNA-seq Data Profiling of the Human Limbus

With the goal of understanding human LESCs and their niche microenvironment, we applied the spatial transcriptomics technique together with scRNA-seq data (downloaded from the database uploaded by De-Quan Li et al. Ocul Surf. 2021;20:20–32. doi: 10.1016/j.jtos.2020.12.004.)⁹ to investigate the human limbus tissues (Fig. 1A). Two healthy human limbal tissues were embedded in the Optimal Cutting Temperature block and sectioned. The tissue sections were fixed, stained with hematoxylin and eosin, and imaged to ensure that the limbus region was included for the following research (Fig. 1B). Subsequently, the sections were subjected to spatial transcriptomics. The number of high-quality spots for space ranger quantitative quality control of each sample is distributed between 444 and 554 and the average number of genes in each spot is distributed between 5057 and 2773 (51 Y for the frontal, and 61 Y for the back). The UMAP of our tissues identified 7 cell clusters (Fig. 1C) and detected the top expressing genes through finding all markers for each cluster and plotting the top 10 genes for each (Fig. 1D). Combined with the histological location, we could see that clusters 1, 2, 4, and 6 were the main components of the limbus region. Based on our results, we observed that keratocan, which is reported to play an important role in the maintenance of corneal homeostasis and function,²⁰ showed the highest expression in cluster 1. SMOC1 is highly expressed in cluster 2 that is close to stroma layer of the limbus, which can encode secreted protein acidic and rich in cysteine–related modular calcium binding 1 and is essential for ocular development in human and mice.^21,22 These results suggest that the components of the stroma layer may play a unique role in regulating the fate of LESCs. 

Figure 1.

 

Overview of spatial transcriptome profiling in the human limbus. (A) Overview of the spatial transcriptomics protocol workflow. (B) Hematoxylin and eosin staining images of the two limbus tissues (51-year-old, left; 61-year-old sample, right). (C) UMAP of two tissues identified seven cell clusters. (D) Heatmap displayed the top 10 genes for each cluster.

To achieve a more precise identification of cell types within each spot, we leveraged previously uploaded scRNA-seq data of human limbal tissue from the database, conducting a comprehensive analysis. A total of 12 clusters were identified by UMAP (Fig. 2A). Different from the original article, we introduced limbal niche cells (mainly including limbal MSCs [LMSCs] and stromal cells) in the analysis²³ and found that cluster 12 highly expressed the niche cell marker gene ALCAM, ICAM1, VCAM1, ITGB1, and THY1^24–26 (Fig. 2B). Immunofluorescence staining of ALCAM was conducted using human limbus tissue sections to further confirm the existence of limbal niche cells (Fig. 2C). Drawing on reported cell representative marker genes, we categorized the 12 clusters into 7 distinct cell types (Fig. 2D). The top three genes for the seven clusters were revealed in Figure 2E. 

Figure 2.

 

Spatial transcriptomics combined with single-cell data reveals heterogeneity and spatial distribution of limbus. (A) UMAP representation of human major limbal classes labeled with 12 cell types. (B) The niche cell marker gene ALCAM, ICAM1, VCAM1, ITGB1, and THY1 are highly expressed in cluster 12. (C) Immunofluorescence staining of ALCAM using human tissues. Scale bar, 50 µm. (D) UMAP of categorizing the 12 clusters into 7 distinct cell types based on cell known marker genes. (E) Heatmap of the top three expressed genes in each cluster. (F) The proportion of each cell type in each spot using RCTD technology. (G) Rank one cell types within each spot in the limbus tissues.

Next, by using the RCTD technology, we estimated the proportion of each cell type in each spot within the tissues and presented it (Fig. 2F). Simultaneously, we illustrated the cell types with the greatest proportion for each spot (Fig. 2G). It could be seen that the surface of the limbus tissue predominantly harbors differentiated cells, conjunctival cells, and progenitor cells. LESCs were localized mainly in the basement membrane, whereas limbal niche cells were situated predominantly within the stromal area. Additionally, melanocytes were found to be dispersed throughout the limbal region. Overall, this comprehensive characterization provides invaluable insights into the cellular distributions of the human limbus, offering a detailed transcriptome map derived from the integration of spatial transcriptomics and scRNA-seq data. 

Different Regions With Unique Gene Expression and Specific Biological Process in the Limbus Area

Previously, a large number of studies had confirmed the x–y–z axis theory of the existence of LESCs in the limbus area.²⁷ To comprehensively unravel the intricacies of the niche microenvironment, we meticulously defined the limbal area into four regions, taking into account the hematoxylin and eosin staining results. These precisely defined regions comprise the inner region, the epithelial layer of the middle region, the stromal layer of the middle region, and the outer region (Fig. 3A). Cell dimension reduction, clustering, and visualization analyses were performed according to specific genes in the region, and four clusters were defined (Fig. 3B). Differential gene expression was analyzed in these clusters, and the top 10 genes were shown by heat map (Fig. 3C). In our results, HES2, HES5, and MUC4 are highly expressed in cluster 4, which close to epithelial layer of middle region. As reported, HES is the downstream transcription genes of the Notch pathway, which is involved in many cellular processes such as cell proliferation and differentiation, as well as the homeostasis of multiple tissues and organs.^2,28,29 MUC4, a component of the glycocalyx, plays an important role in determining the relationship between the cornea and the tear film.^30,31 Taken together, these unique genes provide new insights into the characteristics of the limbus and cross-talk between LESCs and other niche cells. 

Figure 3.

 

Four regions with unique gene expression in various clusters. (A) Four regions of the limbal area according to the hematoxylin and eosin staining results. (B) The cell dimension reduction, clustering, and visualization analysis define four clusters. (C) Heatmap of the top 10 DEGs from each cluster. (D–F) Heatmap for expression of top 15 DEGs in each region.

Subsequently, the DEGs among the four limbal regions on the two tissue sections were analyzed. Remarkably, when contrasted with the other three regions, the epithelial layer of the middle region emerged as the highly expressed genes related to LESCs, such as TP63, KRT5, KRT14, KRT15, GPHA2, and PAX6^32–34 (Figs. 3D–F). This compelling observation strongly suggests that the LESCs predominantly inhabit this specific region, corroborating our understanding of the corneal biology. 

Except for the marker genes previously reported, we found that GPHB5 and PERP were also highly expressed in epithelial layer of the middle region (Figs. 3D–F, Supplementary Fig. S2). As reported, GPHB5 is a secretory protein that has been found to be associated with glucose and lipid metabolism.³⁵ Here, we used immunofluorescence staining to determine whether GPHB5 is expressed in primary LESCs. Surprisingly, the GPHB5 exhibited excellent colocalization with the classic LESC marker, ΔNp63 (Fig. 4A). In addition, we further validated the existence and spatial distribution of GPHB5 in the human limbus tissue by performing multicolor immunofluorescence staining (Fig. 4B), finding that GPHB5 was primarily located in the epithelial layer of the limbus and colocalized with ΔNp63. In addition, to test whether GPHB5 influences the stemness of LESCs, we used silencing RNA to knock down GPHB5 and found the decreased gene expression of stemness marker CK14 and CK15 in si-GPHB5 group, suggesting that GPHB5 may play an important role in the stemness maintenance of LESCs (Fig. 4C). In addition, immunofluorescence staining in LESCs and corneal epithelial cells was also performed. We innovatively found that GPHB5 and other stemness marker such as ΔNp63 and CK14 were highly expressed in LESCs, but barely expressed in corneal epithelial cells (Fig. 4C). Therefore, it is reasonable to speculate that GPHB5 may serve as a novel potential biomarker of LESCs, which has never been reported previously. PERP, a tetraspan membrane protein whose expression is controlled by P63, has an important role in epithelial integrity in many tissues and promotes epithelial cell adhesion.^36,37 By performing multicolor immunofluorescence staining, we innovatively discovered that PERP was highly distributed in the superficial epithelial cells and co-located with corneal epithelial markers CK12 (Supplementary Fig. S3). We inferred that PERP might be a potential candidate for corneal epithelial cell biomarker and the function of them needed to be further verified. Collectively, these DEGs could have guiding significance for understanding the limbal microenvironment and investigating new LESCs and corneal epithelial cells biomarkers. 

Figure 4.

 

Representative immunofluorescence staining images of GPHB5. (A) Immunofluorescent staining co-analysis of GPHB5 with ΔNp63 in LESCs. Scale bar, 100 µm. (B) Immunofluorescent staining of GPHB5 and ΔNp63 in human limbus. Scale bar, 50 µm. (C) The relative gene expression of CK14 and CK15 in the GPHB5 and si-GPHB5 groups. (D) Immunofluorescent staining of GPHB5, ΔNp63 and CK14 in LESCs and corneal epithelial cells. Scale bars, 50 µm. (n ≥ 3, ** P < 0.01, *** P < 0.001.)

Compared with the epithelial layer of the middle region and outer region, RARRES1 and decorin are highly expressed in the stroma layer of the middle region (Fig. 3F, Supplementary Fig. S4A). KRT12 and aldehyde dehydrogenase 3A are highly expressed in the inner region compared with the outer region (Supplementary Fig. S4B). In the outer region, compared with the other three regions, CHI3L1 is highly expressed, which is related to the trabecular meshwork³⁸ (Fig. 3E, Supplementary Figs. S4A and S4B). Taken together, these DEGs might give guidance for understanding the unique development role of these regions. 

Subsequently, GO enrichment (Figs. 5A–C) analysis of DEGs was performed using R (version 4.0.3) based on the hypergeometric distribution. Compared with all three other regions, GO term “Protein Binding” was highly enriched in the epithelial layer of middle region. Protein binding is a crucial process that proteins form bonds with other substances, including other proteins, molecules, and drugs, providing insights to infer the cross-talk of different limbal cells and novel binding sites as targets for therapeutic intervention.³⁹ Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis was also conducted to identify the most over-represented biological pathways in a list of genes (Figs. 5D–F). Based on enrichment score, “ECM–receptor interaction” was the top active in the epithelial layer of the middle region compared with the outer region and the stroma layer of the middle region. Specific interactions between cells and the ECM are mediated by transmembrane molecules, mainly integrins. These interactions directly or indirectly control cellular activities such as adhesion, migration, differentiation, proliferation, and apoptosis. In addition, integrins act as mechanoreceptors and provide a physical link between the ECM and the cytoskeleton that transmits force. Furthermore, some unique biological processes and pathways in other three regions are also demonstrated (Supplementary Figs. S5A–S5E). 

Figure 5.

 

The unique specific biological process, pathways and regulons of four regions in the limbus area. (A–C) Representative GO terms for DEG of each region. (D–F) Kyoto Encyclopedia of Genes and Genomes enrichment of different regions. (G) The top three gene expressions in the four regions. (H, I) Connection specificity index heatmap of the regulon module.

In conclusion, the epithelial layer of the middle region, where LESCs are located, plays an important role in the maintenance of LESCs and communication with other limbal cells to regulate the fate of LESCs. The niche cells are mainly located in the stroma layer of the middle region and are vital in regulating the function of LESCs through various pathways, such as “ECM–receptor interaction.” The inner region is close to the cornea and may be related to the regulation of the proliferation and differentiation of corneal epithelial cells, as well as maintaining corneal transparency. In the outer region, some unique genes that have close connections with the trabecular meshwork and retina are enriched in this region that needed to be further investigated. Taken together, these unique genes and biological processes and pathways provide new inspiration for future studies on the regulation of LESCs, as well as the composition of the limbus. 

Identifying Important Transcriptional Factors and the Regulons of the Four Regions

SCENIC analysis was applied to provide valuable insights into the molecular mechanisms governing the regulation and differentiation of various regions within the limbus. TFs refer to protein molecules that can bind to genes in specific sequences to ensure specific temporal and spatial expression of target genes. Many TFs act as master regulators and select genes that control cell-type determination, developmental patterns, and specific pathways. Here, we predicted changes in the highly active TFs in the isoforms of each block, which were further validated by heatmap expression analysis (Fig. 5G, Supplementary Figs. S6A, S6B). Notably, TP63, a well-studied classical TF and marker gene for corneal development and stem cell proliferation, emerged as a highly active TF enriched in the epithelial basal layer. This finding further confirmed the existence of LESCs in this region. Furthermore, our results brought attention to two additional TFs in this region, namely the GRHL2 and KLF5. GRHL2 assumes an essential TF in inhibiting the epithelial–mesenchymal transition and maintain the homeostasis of human corneal epithelium.⁴⁰ In parallel, KLF5 is essential for the normal maturation and maintenance of the mouse ocular surface and has been proved to regulate a series of genes associated with diverse functions such as cell adhesion, barrier function, and hydration maintenance.⁴¹ These highly expressed TFs may be valuable for future research on the function of LESCs or provide new therapeutic targets for the treatment of LSCD. 

MEF2D and MYCN are highly expressed in the stroma layer of the middle region (Supplementary Figs. S6A and S6B). MEF2D belongs to the MEF2 family, which plays an important role in cell proliferation, differentiation, and death by linking to calcium-dependent signaling pathways.⁴² It is reported that MEF2D is related to the function of adipose-derived MSCs in chronic kidney disease.^42,43 MYCN is one of the most important genes in neuroblastoma research and is essential for retinoblastoma cell growth and tumor formation.⁴⁴ In addition, an early study revealed that MYCN is exclusively expressed in dental mesenchymal cells at E15.5.⁴⁵ Although these two TFs have not been reported to be related directly to the limbus, we speculate that they may also have something to do with the niche cells’ function and thereby regulating the fate of LESCs. 

When it comes to the inner region, interferon regulatory factor 1 (IRF1) is a TF that regulates self-renewal and stress responses, supports the maintenance of hematopoietic stem cells, and plays a central role in the corneal innate immune response by regulating CXCL10 expression.^46,47 KLF has been shown to regulate the growth, proliferation, and differentiation of different cell lineages and is involved in eye development and homeostasis.⁴⁸ In addition, KLF12 was detected to have a specific expression pattern in the cornea.⁴⁹ This further confirms our definition of the inner region, which is close to the cornea. In the outer region, serum response factor is important for limbus angiogenesis during development and the continuing maintenance of the vascular tone.⁵⁰ Their roles in the limbus still need further investigation (Supplementary Figs. S6A and S6B). 

Next, SCENIC software was used to identify modules co-expressed between TFs and regulons to obtain the specific correspondence between the predicted regulon and each cell type (Fig. 5H). The connection specificity index of a regulon was used to represent the association between different regulons. From our results, regulons, such as GRHL2, BCL11A, KLF5, and TP63, are enriched and have a higher connection specificity index in epithelial layer of the middle region (Figs. 5H and 5I), which might imply a critical role in jointly regulating downstream genes and being responsible for cell functions. 

Characterizing Essential Cells in the Limbus

To explore the intricate cell–cell communications within the limbus, we used a specific set of marker genes (Table) to identify four major cell types: LESCs, LMSCs, melanocytes, and immune cells. In addition, immunofluorescence staining and electron microscopy were used to reveal the microstructure and essential cells of the limbus (Supplementary Figs. S7 and S8). This approach allowed us further shedding light on their roles and interactions within the niche microenvironment. 

Table.

 

Specific Marker Genes for Four Essential Cell Types

 

The Cyclone function in the Scran package (version 1.14.3) was used to predict the cell cycle through the expression of “marker gene pairs” (Fig. 6A). Cell-cycle analysis indicated that 20% and 15.4% of LESCs were in the G1 phase in the 51-year-old and 61-year-old samples, respectively. We infer this may be related to age, but a larger sample size is needed for further verification. Interestingly, the immune cells were all in the S phase, which may suggest that they are actively proliferating or replicating their DNA. This could be related to some growth factors inside the limbus; in addition, immune cells might play pivotal roles in regulating the LESCs. 

Figure 6.

 

The cell cycle, gene expression, and intercellular communication analysis of niche cells. (A) The cell-cycle stages for the four limbal cells were predicted by subpopulation using the Cyclone function in Scran (version 1.14.3) package. (B) Expression patterns of top three genes of each cell. (C) Spatial localization of the four major cell types (the repeat represents three or more cell replicates). (D) Heatmap demonstrating the total number of interactions between cell types. (E) Bubble plot of significant interactions between LESCs and other cells. (F) Immunofluorescent staining of CD44 and COL6A2 in human limbus. Scale bar, 50 µm. (G) Chord diagram of COLLAGEN, PECAM1 and CCL pathway. (H) Information flow of each signaling pathway.

Indeed, SCENIC analysis was used to evaluate the regulatory networks within these four cell types, particularly focusing on the activity of TFs and their related regulons (Fig. 6B). The IRF family, including IRF1, IRF6, and IRF7, are important regulators of stress hematopoiesis and normal hematopoiesis.⁵¹ It is notable that the expression of IRF1 is the highest in LESCs. It can regulate the self-renewal and anti-apoptosis ability of stem cells and support the maintenance of hematopoietic stem cells.⁵² IRF7 is highly expressed in LMSCs, and it is reported that CXCR4 is the downstream target of IRF7.⁵³ Interestingly, LESCs interact with MSCs via SDF-1/CXCR4, an axis mainly expressed in the human limbus rather than the cornea, with SDF-1 preferentially expressed by limbal epithelial cells and CXCR4 expressed by MSCs.⁵⁴ Therefore, we speculate that IRF7 may play an important role in the communication between LESCs and MSCs, which requires further investigation. In melanocytes, the expression of IRF6 was elevated. The methylation levels of 5′ IRF6 CpG island are potentially associated with the sensitivity of melanoma cells to interferon. IRF6 also plays a tumor suppressor role in keratinocytes by mediating Notch differentiation.⁵⁵ 

Except for the above, the expression of activating TF 2 (ATF2), ATF3, and ARID5B was increased in immune cells. ATF2 and ATF3 belong to the ATFs, which play an important role in cell proliferation, apoptosis, differentiation and inflammatory related pathological processes.⁵⁶ The AT-rich interaction domain (ARID) protein family, especially ARID5B, and are widely expressed in adult organs, including the lungs, small intestine, kidneys, muscles, heart, and brain.⁵⁷ Importantly, ARID5B affects the development and function of mouse B cells, while playing an important role in regulating NK cell metabolism.^58,59 Taken together, these new TFs provide a better understanding of the function and characteristics of LESCs and lay a foundation for exploring the pathogenesis of LSCD. 

At the same time, we also showed the spatial localization of the four major cell types based on a specific set of marker genes which mentioned above (Fig. 6C, Supplementary Figs. S9A–S9D) to better explore their cross-talk at the spatial level. We observed that there were a large number of overlapping spots of different cells in the epithelial layer of the middle region, especially between LESCs and LMSCs. In the outer region, the colocalization of LESCs with other cells was reduced, and LMSCs were at deeper levels. 

Niche Regulation of LESCs Revealed by Cell–Cell Communications

To reveal the mechanisms of niche regulation at the human limbus, we attempted to explore intercellular communication among the four essential cells mentioned above. 

On the basis of the ligand–receptor analysis, we determined and counted the number of interacting links between the four cell types (Fig. 6D). It is noteworthy that, when LESCs act as receptor cells, LMSCs have the highest number of ligands associated with LESCs, with melanocytes following closely behind LMSCs. This finding proves that, in the limbal microenvironment, LMSCs have the most important significance in regulating the fate of LESCs. 

Next, the CellChat R package (version 1.1.3) was used to analyze cell communication between the outer region and the middle region (including epithelial layer of the middle region and stromal layer of the middle region). It is noteworthy that COL6A2/CD44 is widely active in immune cells and LESCs, melanocytes and LESCs, LESCs and LESCs, and LMSCs and LESCs through ECM–receptor interaction by targeting COLLAGEN pathway (Fig. 6E). Multicolor immunofluorescence staining was conducted to validate the activation of COL6A2/CD44 between limbal stem cells and other cells (Fig. 6F). In addition, PECAM1/PECAM1 and CCL14/ACKR1 are active in immune cells and LESCs, LESCs and LESCs, and LMSCs and LESCs through cell–cell contact by targeting the PECAM1 pathway and through sSecreted signaling by targeting CCL pathway, respectively (Figs. 6E and 6G). 

In addition, LAMB3/CD44 and MIF/CD74-CD44 are active in melanocytes and LESCs through ECM–receptor and secreted signaling, respectively. FN1/CD44 is active in LMSC and LESCs through ECM–receptor interaction. COL4A2/CD44 is active in immune cells and LESCs through the ECM–receptor interaction (Fig. 6E). Moreover, except for the pathways mentioned, we observed that the MIF and LAMININ pathways are active in the middle region and the NOTCH, CXCL, and LAMININ pathways are active in the outer region (Fig. 6H). Taken together, these results suggest a complex interaction between LESCs and their niche cells, and these ligands and receptors, known and unknown, as well as these unique signaling pathways, may provide new opportunities to explore the mechanism of regulating LESCs. 

Discussion

The limbus maintains the homeostasis and angiogenic privilege of cornea, but little is known about how the niche microenvironment regulates the LESCs. The deep exploration of the limbus at the spatial level can provide new insights for understanding LESCs and surrounding components. In this study, we generated a spatial atlas of the human limbus and located almost all major cell types for the first time. By using various bioinformatics analyses, we demonstrate the functional characteristics of different regions within the limbus. Furthermore, we elucidate the intricate interactions between the four major cell populations inside the limbus. These findings not only deepen our understanding of limbal biology, but also lay the foundation for future investigations into LESCs and the development of therapies for limbal stem cell deficiency. 

In this study, we used scRNA-seq data from human LESCs uploaded by Chen et al. Different from the original article, we found out that the cluster 12, as described by Li et al., exhibited high expression levels of several marker genes associated with limbal niche cells, including ALCAM, ICAM1, VCAM1, ITGB1, and THY1. Based on this observation, we hypothesized that cluster 12 represents a population of limbal niche cells and proceeded to map it onto our spatial transcriptomics results for further analysis. Despite cluster 12 comprising only 22 cells according to Chen's research, and being broadly distributed across spatial spots, we attribute this discrepancy to the fact that the limbal epithelial tissue samples analyzed in Chen's study may have contained a relatively limited number of niche cells. Consequently, variations in gene expression profiles could arise due to the scarcity of these niche cells within the sampled tissue. By integrating the single-cell data, we effectively mapped the spatial distribution of cells within the limbal tissue. This comprehensive approach not only provided us with a detailed spatial representation of cellular composition, but also facilitated a deeper understanding of the limbal microenvironment. 

At present, several research have confirmed the x–y–z axis theory of LESCs, but the specific functions of different regions of the limbus are still not well-studied. Therefore, four regions of the limbus were divided innovatively based on histological structure, and the DEGs, as well as the biological process and signal pathways, were analyzed thoroughly. In our study, we found that GPHB5 and PERP are highly expressed in the epithelial layer of the middle region based on the results of immunofluorescent staining. Compared with the epithelial layer of the middle region and outer region, RARRES1 and decorin are highly expressed in the stroma layer of the middle region. RARRES1 is not only an important regulator of autophagy, but also an inhibitor of retinoic acid–regulated carboxypeptidase, which is related to aging, metabolism, and stem cell differentiation.^60,61 Whether it has any effect on LESCs needs further investigation. Decorin is a small leucine-rich proteoglycan that mainly exists in the ECM and is involved in the regulation of corneal wound healing.⁶² Compared with the outer region, KRT12 and aldehyde dehydrogenase 3A are highly expressed in inner region. Aldehyde dehydrogenase 3A1 is one of the most abundant proteins expressed in the corneal epithelium and can regulate the proliferation and differentiation of corneal epithelial cells, as well as maintain corneal transparency and protect intraocular tissues from oxidative damage caused by ultraviolet radiation.^63,64 Collectively, these DEGs might give guidance for understanding the unique development role of these regions. 

At the same time, “ECM–receptor interaction” was both activated in the epithelial layer of middle region and the stromal layer of middle region. As reported, paracrine factors and their receptors, cell–cell contacts, cell–matrix contacts, and mechanical transduction are important for LESC self-renewal and fate determination. Inner region is close to cornea and may be related to the regulation of the proliferation and differentiation of corneal epithelial cells, as well as maintaining corneal transparency. These unique findings lay a foundation for exploring the characteristic of the limbus and provide a new perspective for figure out the regulation of LESCs at spatial level. However, some genes have not been reported or researched in depth, further studies are still needed. 

The corneal limbus hosts a consortium of MSCs, melanocytes, and immune cells, collectively exerting crucial regulatory influences over LESCs’ fate. Furthermore, LESCs themselves are integral contributors to this regulatory microenvironment. Indeed, beyond the secretion of multiple growth factors, the intricate network of cell–cell contacts among these cellular constituents also plays a critical role in the regulation of LESCs. Presently, spatial transcriptomics emerges as an invaluable tool for elucidating the intricate signaling among these cells, thereby for the first time, we applied this technology for exploring the niche microenvironment. In this study, we described four major cells based on some marker genes: LESCs (TSPAN7, SOX17, SELE, ECSCR, RAMP3, and IGFBP4), LMSCs (ENG, ALCAM, ICAM1, and VCAM1), melanocytes (TYR, TYRP1, MLANA, and MC1R), and immune cells (CD4, CD44, FOXP3, or NCR1), and emphasized the intercellular communication between four major cell types to further explore the cross-talk in the limbal microenvironment. The reasons why we have not selected more marker genes are that each spot is a combination of multiple cells, and too many genes will lead to a higher proportion of false negatives. 

In the analysis of CellChat, we combined epithelial layer of the middle region and stromal layer of the middle region with the middle region, because there were too few spots in each to analyze. At the same time, there are almost no positive spots in the inner region, so we exclude the inner region and mainly analyzed the CellChat between the middle region and the outer region. The main difference between the middle region and the outer region is that a large number of LESCs enriched in the middle region, so the differences in subsequent analyses are mainly due to the variation of LESCs. Based on our results, we found that LMSCs have the highest number of ligands associated with LESCs, with melanocytes following closely behind LMSCs when LESCs act as receptor cells. This finding proves that, in the limbal microenvironment, LMSCs play a pivotal role in the regulation of LESC fate. In addition, we observed that LMSCs were located mainly in the stromal layer of the middle region and the outer region based on immunofluorescence staining. As in previous reports, LMSCs interact with LESCs through a variety of pathways, including IL-6/STAT3, SDF-1/CXCR4, BMP/Wnt, vimentin, and aquaporin-1.^65,66 In our study, we innovatively found that ligand receptors, such as COL6A2/CD44 and FN1/CD44, were active in LMSC and LESCs through ECM–receptor interaction and PECAM1/PECAM1 is active through cell–cell contact, which may be beneficial for the subsequent study of the interaction between LMSC and LESCs. It is noteworthy that COL6A2/CD44 is also widely active in immune cells and LESCs, melanocytes and LESCs, and LESCs and LESCs through ECM–receptor interaction. Multicolor immunofluorescence staining was conducted to validate the activation of COL6A2/CD44 between limbal stem cells and other cells. Melanocytes, which are located in the basal layer of limbal epithelium, are thought to protect LESCs against UV damage and maintain the phenotype of LESCs. Except for the ligands and receptors mentioned above, it was discovered that LAMB3/CD44 and MIF/CD74-CD44 are also active in Melanocyte and LESCs through ECM-receptor interaction and secreted signaling, respectively. Immune cells can promote stem cell proliferation in various tissues, but few similar observations have been reported in the cornea and there are fewer studies to figure out the interactions between immune cells and LESCs. In our results, we innovatively found that COL4A2/CD44 and PECAM1/PECAM1are active in immune cells and LESCs through ECM–receptor interaction and cell–cell contact respectively, which lay a foundation for future study in immune cells in the limbus. At the same time, PECAM1/PECAM1 and COL4A1/CD44 are active in LESCs and LESCs through cell–cell contact and ECM–receptor interaction, respectively. Taken together, these results suggest a complex interaction between LESCs and their niche cells, and these ligands and receptors, known and unknown, as well as these unique signaling pathways, may provide new opportunities and data support to explore the mechanism of regulating LESCs. 

At the same time, in the 51-year-old sample and 61-year-old sample, we found that the proportion of LESCs in G1 phase gradually decreased, so we hypothesized that the cycle of LESCs, as well as their state, was related to age. However, it still needs more sample size support and follow-up research. 

There is no denying that this study has certain limitations owing to the small sample size, narrow age range, and the lack of further verification of the function of some unique genes. In addition, the donor cornea will affect the results of the study, such as age, cause of death, time of absence, location, and size of the limbal collection. In addition, the 10× Genomics Visium platform used in this study suffered from some technical limitations and the accuracy of this platform needed further improvement. However, the biological information provided by this technology can undoubtedly help researchers better understand the spatial structure of the limbus and the signal transduction therein. 

In summary, our analysis reveals the spatial mapping of the limbal region and the deeper interactions between various niche cells, and forms a powerful translational resource that provides important insights into the limbal niche of LESCs and new strategy for treating LSCD. At the same time, this dataset can be used as a baseline for the limbal niche in healthy humans for the assessment of cellular composition and transcriptomic analysis, providing a searchable resource for the scientific and clinical community. 

Acknowledgments

Supported by the National Natural Science Foundation of China (82201136, 82271041, 82070919, and 82471040), the Program of Shanghai Academic/Technology Research Leader (22XD1401800), the Biomaterials and Regenerative Medicine Institute Cooperative Research Project, Shanghai Jiao Tong University School of Medicine (2022LHA06), the Disciplinary Crossing Cultivation Program of Shanghai Jiao Tong University (YG2022QN055), Shanghai Key Clinical Specialty, and Shanghai Eye Disease Research Center (2022ZZ01003). 

Author Contributions: S.L. and L.C. drafted the manuscript and created all the figures. H.S. and F.F. assisted in vitro experiments. S.Z. and J.C. examined the data. S.L., H.S., F.F., C.S., Y.F. and L.C. discussed the concepts of the manuscript. All authors read and approved the final manuscript. 

Disclosure: S. Li, None; Hao Sun, None; F. Fang, None; S. Zhang, None; J. Chen, None; C. Shao, None; Y. Fu, None; Liangbo Chen, None