You are viewing the site in preview mode

Skip to main content
Download PDF

Spatiotemporal single-cell analysis elucidates the cellular and molecular dynamics of human cornea aging

Genome Medicine volume 17, Article number: 56 (2025) Cite this article

Abstract

Background

The human cornea is a transparent and uniquely ordered optical-biological system. Precise coordination of its cellular mechanisms is essential to maintain its transparency and functionality. However, the spatial, cellular and molecular architecture of the human cornea and its intercellular interactions during aging have not been elucidated.

Methods

We performed single-cell RNA sequencing (scRNA-seq) and single-cell SpaTial Enhanced REsolution Omics-sequencing (scStereo-seq) analysis in corneal tissue from eight eyes of donors aged 33-88 years to elucidate the spatiotemporal cellular and molecular dynamics of human cornea aging. Immunofluorescence staining and Western blotting were performed to validate the findings.

Results

Spatiotemporal single-cell analysis revealed the complex cellular landscape, spatial organization and intercellular communication within the human cornea. The subpopulations of major cell types of the cornea were elucidated with precise spatial positions. In particular, we identified novel subpopulations, mapped the spatial positioning of limbal stem cells within the limbal niche, and delineated the interactions between major cell types. We observed that three basal cell subsets migrate centripetally from the peripheral to the central cornea with age, suggesting the “spatiotemporal centripetal pattern” as a novel paradigm for the age-related migration of corneal epithelial cells. Furthermore, we elucidated the age-related, region-specific molecular and functional characteristics of the corneal endothelium, demonstrating differential metabolic capacities and functional properties between the peripheral and central regions.

Conclusions

As the first comprehensive spatiotemporal atlas, our work provides a valuable resource for understanding tissue homeostasis in the human cornea and advances research on corneal pathology, transplantation, senescence and regenerative medicine in the context of corneal aging.

Background

The cornea is a unique transparent avascular tissue that maintains visual clarity and integrity [1,2,3]. Corneal homeostasis, which is critical for transparency, requires a dynamic balance of cellular processes that can be disrupted by aging [4]. The human cornea has a stratified architecture comprising three cellular layers: the corneal epithelium, the stroma, and the corneal endothelium. The outermost epithelium is continuously renewed by limbal stem cells (LSCs) within a specialized niche, with a full cycle of cellular regeneration typically completed in approximately 7–14 days. Keratocytes, the mitotically quiescent cells of the corneal stroma, can transition into a myofibroblastic state upon injury, contributing to the formation of fibrotic scars [5]. The endothelium, a monolayer of flattened hexagonal cells, lacks regenerative capacity and declines with age, which may affect its function [6]. Given the complexity and delicate regulation required to maintain the dynamic homeostasis of the cornea, an understanding of its cellular composition as well as its spatial and molecular characteristics is essential to gain insight into corneal health and disease.

Current knowledge of the cellular and molecular architecture of the human cornea, particularly how it changes with advancing age, remains incomplete. Despite advances in single-cell RNA sequencing (scRNA-seq) that have revealed cellular heterogeneity in human limbal tissue [7,8,9], traditional scRNA-seq analyses, while informative, often lack the spatial context necessary to fully understand intercellular interactions and regional differences within the cornea [10, 11]. Furthermore, the paucity of whole human corneas from different age groups and the consequent reliance on animal models or limbal tissues have limited our understanding. This knowledge gap hinders advancements in corneal therapy, transplantation, and regenerative medicine. The development of single-cell SpaTial Enhanced REsolution Omics-sequencing (scStereo-seq) now allows high-resolution, nanoscale analysis of the spatial and functional organization of cells in tissues [12,13,14].

This study systematically characterized the cellular composition and spatial organization of major corneal cell types in whole human corneas and adjacent conjunctiva from donors aged 33–88 years at single-cell resolution. By combining scStereo-seq, scRNA-seq, Western blot and immunofluorescence, we provide a detailed spatiotemporal transcriptomic and cellular atlas of the whole human cornea. This atlas is a valuable resource for advancing research in corneal pathology, transplantation, and regenerative medicine, and provides new insights into the spatiotemporal hallmarks of human corneal aging.

Methods

Ethical declaration

This study was conducted in accordance with the principles of the Declaration of Helsinki. The use of corneal tissues was approved by the Ethical Review Committee of the Affiliated Eye Hospital of Wenzhou Medical University (2022-159-K-124-01). Written informed consent was obtained from the legal representatives of the donors before tissue collection. The mouse experiments were approved by the Laboratory Animal Ethics Committee of Wenzhou Medical University (YSG24032101).

Human specimens

The corneal tissues used in this study were collected from four male donors aged 33, 50, 68, and 88 years (n = 8 corneas), each including a 3-mm rim of sclera and conjunctiva. These samples were surgically excised from donor eyes obtained from the Eye Bank at Wenzhou Medical University. All protocols were in accordance with the Interim Measures for the Administration of Human Genetic Resources, which are administered by the Ministry of Science and Technology of China. The donor criteria included absence of ocular pathology or systemic diseases that could affect the cornea, and death from causes unrelated to eye disease. The post-mortem time to tissue collection was less than 24 h. Prior to inclusion in this study, all donor corneas underwent morphological examination and fluorescein sodium staining to confirm corneal epithelial integrity. The human corneal samples were processed for scRNA-seq and scStereo-seq analyses, histological examination by hematoxylin and eosin (H&E) staining, and immunofluorescence staining for multiple marker proteins.

Animal samples

Male C57BL/6 mice aged 8 weeks and 14 months (6 for each group) were obtained from and housed in the specific-pathogen-free facility at Wenzhou Medical University under standard conditions with controlled temperature (22 ± 2°C), humidity (50 ± 10%), and a 12-h light/dark cycle. Mice were provided with standard laboratory chow and water ad libitum. The cages were changed twice weekly, and animal health was monitored daily. Mice were euthanized by cervical dislocation following inhalation of 2% isoflurane for initial anesthesia, in accordance with the guidelines for the humane treatment of laboratory animals. Immediately after euthanasia, corneal tissues were collected for subsequent analyses. For Western blotting, corneas from young and aged mice (3 for each group) were dissected into central and peripheral regions to examine neurotrophic receptor tyrosine kinase 2 (NTRK2) expression. For whole-mount immunofluorescence staining, intact corneas from young mice (n = 3) were carefully excised to examine the expression patterns of NDUFV2 and succinate dehydrogenase B (SDHB) in the corneal endothelium.

scRNA-seq sample preparation and sequencing

Whole corneas with 3 mm white sclera and conjunctiva were dissected from donor eyes. The tissue was transferred directly to a clean RNase-free culture dish and separated into fragments. The fragments were then dissociated into individual cells via a dissociation solution in a water bath at 37 °C with agitation at 100 rpm for 20 min. The resulting cell suspension was filtered through a 20-μm cell strainer and then centrifuged to collect the cells. Erythrocytes were removed via lysis buffer, followed by resuspension of the cell pellet and removal of nonviable cells. The final cell pellet was resuspended in 50 μl of 1 × PBS with 0.04% bovine serum albumin, and cell viability was confirmed by trypan blue exclusion. A viability over 85% was considered acceptable, and the cell concentration was adjusted to 700–1200 cells/μl. Single-cell suspensions were then obtained via the 10X Genomics Chromium Single-Cell 3′ Kit according to the manufacturer’s instructions, and the captured cells were subsequently amplified and processed for library construction according to a standard protocol. Libraries were sequenced on an Illumina NovaSeq 6000 platform (Illumina, San Diego) at 150 bp paired-end read lengths by Shanghai Personal Biotechnology (Shanghai, China).

scStereo-seq sample preparation and sequencing

For scStereo-seq, donor corneas with 3 mm sclera and conjunctiva were excised from the isolated eyeball. The sutured tissues were embedded at the optimal cutting temperature (OCT), precooled with liquid nitrogen and isopentane, and stored at − 80 °C prior to the experiment. The tissues were sectioned at a thickness of 10 μm via a Leica CM1950 cryostat, attached to the scStereo-seq capture chip, and warmed for 3 min at 37 °C in a slide warmer. The chip was then incubated at − 20 °C in methanol for 30 min, stained with nucleic acid dye (Thermo Fisher, Q10212) and imaged (Ti- 7 Nikon Eclipse microscope). The cDNA products were extracted and purified from the scStereo-seq chips after permeabilization and reverse transcription. Indexed cDNA libraries were then prepared according to the manufacturer’s instructions and sequenced on an MGI DNBSEQ-Tx sequencer (read 1: 35 bp, read 2: 100 bp).

H&E staining

Corneal tissues were embedded in OCT compound, and frozen sections were cut into 10 μm slices for H&E staining. The sections were first fixed in 4% paraformaldehyde (PFA) for 10 min and then washed twice for 2 min with distilled water. The sections were then stained with hematoxylin for 5–10 min, immersed in tap water for 10 min, and washed again with distilled water. Eosin counterstaining was then performed for 1 min. The corneal sections were then dehydrated sequentially in 70%, 80%, 90%, and 100% ethanol for 10 s each and washed with xylene twice for 5 min each.

Immunofluorescence staining

The slides were fixed in 4% PFA for 15 min at room temperature and washed three times in PBS. Fixed sections were permeabilized in 0.4% Triton X- 100 in PBS for 20 min and then blocked for 2 h with 10% normal goat serum at room temperature. The tissue sections were incubated with primary antibodies overnight at 4 °C, followed by fluorescence labeling with secondary antibodies and immunoimaging via a confocal microscope. The following antibodies were used for immunofluorescence: anti-KRT15 (Abcam, ab52816, 1:200), anti-IFITM3 (Invitrogen, PA5 - 11,274, 1:200), anti-HIF1A (Proteintech, 66730–1-Ig, 1:200), anti-UBE2S (Invitrogen, PA5 - 142029, 1:200), anti-CLDN4 (Invitrogen, 329494, 1:200), anti-SLC6A6 (Abcam, ab236898, 1:200), anti-KRT15 (Santa Cruz, sc70912, 1:200), anti-NTRK2 (Santa Cruz, sc-136991, 1:200), anti-CCL14 (Proteintech, 14216–1-AP, 1:200), anti-LY6D (Proteintech, 17361–1-AP, 1:200), anti-CA3 (Proteintech, 15197–1-AP, 1:200), anti-sodium potassium ATPase (Abcam, ab283318, 1:200), anti-goat IgG (Alexa Fluo 647, Invitrogen, 1:400), anti-rabbit IgG (Alexa Fluo 555, Invitrogen, 1:400), anti-rabbit IgG (Alexa Fluo 594, Invitrogen, 1:400), anti-mouse IgG (Alexa Fluo 594, Invitrogen, 1:400), anti-mouse IgG (Alexa Fluo 488, Invitrogen, 1:400), anti-mouse IgG2b (CoraLiteb Plus 555, Proteintech, 1:500), anti-mouse IgG1 (Alexa Fluo 647, Proteintech, 1:500), and anti-mouse IgG2a (CoraLiteb Plus 488, Proteintech, 1:500).

Western blotting

Corneal tissues were lysed in radioimmunoprecipitation assay buffer (Solarbio, Beijing, China) and protein concentrations were measured using a bicinchoninic acid protein assay kit (Beyotime, Shanghai, China). Equal protein amounts were separated by sodium dodecyl-sulfate polyacrylamide gel electrophoresis (GE Healthcare Life Sciences, USA) and transferred to Polyvinylidene Fluoride membranes. Membranes were blocked with 5% milk for 2 h at room temperature, then incubated with primary antibody against NTRK2 (Santa Cruz, sc-136991) overnight at 4 °C. After washing, membranes were incubated with HRP-conjugated secondary antibodies for 2 h. Protein bands were visualized using chemiluminescence reagents (Beyotime, China) and imaged using a GE Amersham Imager AI680. Band intensities were quantified using ImageJ software and normalized to loading controls.

scRNA-seq data quality control and preprocessing

For the scRNA-seq data, the raw reads were processed using Cell Ranger (v7.0) with default parameters and aligned to the GRCh38 genome reference and annotation. The gene expression matrix for each cell was then generated and imported into a Seurat object via the Seurat R package (v4.0.0) for downstream analyses [15]. For each sample, only cells with less than 10% mitochondrial counts and more than 200 detected genes were retained, and genes expressed in more than three cells were selected for further analysis. Potential doublets were detected and removed via the Scrublet pipeline (v0.2.3) package, with an expected_doublet_rate of 0.06 [16]. After strict quality control, we ultimately obtained 25,604 cells. The SCTransform() function was used to normalize the variance in sequencing depth and to regress the mitochondrial percentage and cell cycle scores. The datasets from different libraries were integrated via the RunHarmony() function for batch effect correction with Harmony (v0.1.1) [17]. Dimensionality reduction was performed via RunPCA() and the top 30 PCs were used for visualization via uniform manifold approximation and projection (UMAP). Cell clustering was conducted via the FindClusters() function with a resolution of 0.1. Active LSCs were clustered into three subsets at a resolution of 0.6, and stromal cells were reclustered at a resolution of 1.5. Cell cycle parameters (G1, S and G2/M phases) were estimated via the Seurat function CellCycleScoring().

Identification of differentially expressed genes

To define cluster-specific markers, we performed differential gene expression analysis between groups or clusters to identify differentially expressed genes (DEGs) based on normalized gene expression values via the Wilcoxon rank-sum test of the FindAllMarkers() function. Only genes with log fold changes > 0.25 and adjusted p values < 0.05 were considered significant for further analysis.

Enrichment analysis

Gene ontology (GO) enrichment analysis for specific gene sets was performed via the Database for Annotation, Visualization, and Integrated Discovery (DAVID) webserver [18] and the gene set enrichment analysis (GSEA) R package with default parameters (v1.2).

Fuzzy cluster analysis

Time-series gene sets were detected via soft clustering, which was performed with the Mfuzz R package [19]. For specific cell subpopulations, we first calculated the average expression values of all genes across different age groups within this subpopulation. Low-expression genes (with an average expression of less than 0.1 in the subpopulation) and genes with low variability (with a standard deviation of less than 0.1 in average expression across age groups) were removed. Genes with the highest average expression in either the 33 years or 88 years were retained for subsequent filtering of aging-related gene sets. The remaining genes were clustered using the mfuzz() function. The optimal number of clusters was determined based on the Hubert statistics values. The gene sets from clustering were ranked according to their membership values, and the top 100 genes were subjected to GO enrichment analysis.

Identification of the cellular neighborhood

To construct cellular neighborhoods (CNs), we used a window-based capture approach, which involves identifying the n nearest neighboring cells surrounding a target cell (n = 10 in our work), as outlined in a previous study [20]. Each window is represented as a frequency vector that captures the composition of the X closest cell types relative to the reference cell. Each cell was then assigned to a CN based on its corresponding window vector.

Pseudotime trajectory analysis

Pseudotime trajectory analysis was performed via the Monocle2 R package (v2.24.1) with default parameters [21]. The trajectories were visualized in a two-dimensional space via dimensionality reduction via the DDRTree reduction method in the reduceDimension function(). Finally, the cells were ordered according to pseudotime via the order_cells function(), and the trajectory was characterized via the plot_cell_trajectory() function.

Transcription factor analysis

The activated transcription factors (TFs) of the cell types were identified via pySCENIC (v0.12.0) [22]. The raw count expression matrix in Seurat was used as input data [23]. GRNBoost was used to infer potential TF targets via default parameters [24]. We identified TFs that were significantly differentially expressed in the four LSC subsets (log fold change > 0.58 and adjusted p-value < 0.05) as their specific TFs and screened target genes with importance greater than 0.1 to construct the TF network via Cytoscape (v3.10.2).

Cell‒cell communication analysis

Cell‒cell communication between cell groups of interest was inferred via CellChat (v1.1.1) with an existing ligand-receptor interaction database [25]. In brief, CellChat identified ligands or receptors that were significantly overexpressed in a cell group. The strength of communication between two interacting subsets was assessed by quantifying the average expression level of ligands in one subset and receptors in the other. The inferred cell‒cell communication network was visualized via the chordDiagram function in the circlize R package (v0.4.13).

Spatial scStereo-seq data processing

Paired-end raw sequencing fastq reads were generated via the MGI DNBSEQ-TX sequencer. Raw in situ spatial transcriptomics data were processed and unsupervised clustered as previously described [13]. The final gene expression profile matrices were obtained via SAW, a pipeline available on GitHub (https://github.com/BGIResearch/SAW) [26]. The sequences with invalid CID sequences (1 base mismatch allowed) and low-quality unique molecular identifier (UMI) sequences (quality score < 10) were removed and then mapped to the reference genome of Homo_sapiens (hg38) via STAR [27]. Transcripts captured in bin 50 X 50 DNA nanoballs (bins) were merged into bin 50 and transferred to Seurat for downstream analyses. The SCTransform function was used to normalize and scale the raw matrix of the scStereo-seq data. We separately computed the top 1500 highly variable genes (HVGs) for each sample, merged these genes as the HVGs for the integrated scStereo-seq data, and used them as inputs for principal component analysis for dimensionality reduction. Finally, we performed unsupervised spot clustering via the FindCluster function to divide spatially distinct regions with a resolution parameter of 0.2.

Cell type deconvolution in the Stereo-seq data

Spatial deconvolution of cell types was performed from matched scRNA-seq data via the robust cell type decomposition (RCTD) algorithm in the spacexr R package (v2.2.1) [28]. For each slice in the scStereo-seq data, spatial coordinates and raw count matrices were used to construct the spatial RNA object in the RCTD. The raw count matrix and cell identity of the scRNA-seq data were used as a reference. The reference and spatial RNA objects were then input into the run. RCTD() function with the doublet mode is used to calculate the weights of each cell subset. Finally, the cell subset with the maximum weight quantified by RCTD was used to annotate each spot.

Gene signature score analysis

Gene signatures related to stem cell differentiation, epithelial cell differentiation, oxidative phosphorylation, ion transport, and ATP biosynthesis processes were collected from MSigDB (https://www.gsea-msigdb.org/gsea/msigdb/) [29]. Collagen, extracellular matrix (ECM) proteoglycans, and elastic fiber-related genes were obtained from Maiti et al. [30]. The gene signature score was calculated for each cell via the AddModuleScore() function in Seurat (v4.0.0).

Construction of the thrombospondin- 1-associated gene module

To construct the thrombospondin-1 (THBS1)-associated gene module, we calculated the Pearson correlation coefficient for each gene relative to THBS1 within the NTRK2+ basal cells and identified 43 genes with correlation coefficients greater than 0.2 as highly correlated with THBS1. These genes were then entered into the STRING database to construct a protein–protein interaction (PPI) network. Finally, we obtained a connected PPI network consisting of 26 genes, with THBS1 as the hub gene, which served as the THBS1-associated gene module.

Mouse corneal whole mounts

To analyze the expression profiles of specific proteins in the corneal endothelium, we performed whole-mount immunofluorescence staining on excised corneas from C57/B6 mice. The corneal tissues were fixed in 4% PFA for 15 min at room temperature and then washed three times with PBS. Next, the tissues were permeabilized in 0.4% Triton X- 100 in PBS for 20 min and blocked with 10% normal goat serum at room temperature for 2 h. We incubated the corneas overnight at 4 °C with the following primary antibodies: anti-NDUFV2 (I) (Proteintech, 15,301–1-AP) and anti-SDHB (II) (Proteintech, 10,620–1-AP). After washing, the corneas were incubated with the corresponding secondary antibodies conjugated to Alexa Fluor dyes (Invitrogen) or CoraLiteb Plus dyes (Proteintech) for 2 h at room temperature. Finally, we obtained Z-stack images of the corneal whole mounts via a confocal microscope (LSM 880 META; Carl Zeiss GmbH, Hamburg, Germany).

Statistical analysis

All the statistical analyses were performed in R (v4.3.1). The Wilcoxon rank-sum test was used to compare two groups. Spearman’s rank correlation coefficient was used for correlation analysis. A two-tailed p-value < 0.05 was considered statistically significant.

Results

Spatiotemporal transcriptomic and cell type atlas of the whole human cornea

To elucidate the cellular diversification and spatiotemporal developmental dynamics of the human cornea, we performed scRNA-seq and spatial transcriptomic analyses of corneal tissue from four healthy male donors via 10X Genomics and scStereo-seq technology (Fig. 1a). We obtained a total of 25,604 cells from all samples following stringent quality control and filtering (Additional file 1: Fig. S1a). Unsupervised clustering revealed 19 distinct cell clusters after batch effect correction (Fig. 1b and Additional file 1: Fig. S1b). These clusters were annotated as ten major cell types based on canonical marker gene expression: conjunctival epithelial cells, corneal epithelial cells, limbal cells, stromal cells, smooth muscle cells, vascular endothelial cells, lymphatic endothelial cells, melanocytes, T cells and Langerhans cells (Fig. 1c). The accuracy of the cell type annotations was confirmed by the identification of upregulated genes and biological features specific to each cell type (Additional file 1: Fig. S1c and Additional file 2: Table S1). Corneal epithelial and limbal cells were consistently identified as the predominant cell types in all samples (Fig. 1d).

Fig. 1
figure 1

Spatiotemporal transcriptomic and cell type atlas of the human cornea. a Schematic illustration of the human cornea used for scRNA-seq (10X Genomics), spatial transcriptomics (scStereo-seq) and immunofluorescence. b UMAP visualization of the scRNA-seq data of the 19 primary cell clusters. c Dot plot illustrating the expression of marker genes for each cell cluster. The color indicates the normalized expression level, while the dot size indicates the proportion of expressing cells. d Stacked plot showing the proportion of cell types in each sample. Immune cells include T cells and Langerhans cells. e UMAP visualization of spatial transcriptomics data showing the seven primary clusters. f Spatial visualization of the seven primary clusters in each slide (top), H&E staining of each human cornea (bottom). g Dot plot illustrating the expression of marker genes for each anatomical region. The color indicates the normalized expression level, while the dot size indicates the proportion of expressing cells. h Schematic illustration of cell type mapping from scRNA-seq data to spatial transcriptomics data (left). Spatial visualization of all the cell clusters in each slide (right)

Full size image

To further elucidate the spatiotemporal developmental dynamics of the human cornea at the single-cell level, we performed spatial transcriptomic analyses via scStereo-seq technology (Fig. 1a). We chose a bin 50 configuration (50 × 50 bins, 25 μm diameter per spot, representing a single predominant cell type) to achieve a balance between spatial resolution and gene coverage [13], resulting in an average of 1123 unique molecular identifiers (UMIs) and 505 genes per spot. A total of 132,520 spots were obtained (Additional file 1: Fig. S1 d). Each spatial transcriptomic sample was highly correlated with its corresponding scRNA-seq sample in terms of overall gene expression (Additional file 1: Fig. S1e). Unsupervised clustering categorized all the spots into seven clusters, each corresponding to a specific anatomical region: the ciliary body, conjunctival epithelium, conjunctival crypt, corneal epithelium, corneal endothelium, sclera, or corneal stroma (Fig. 1e). These clusters correlated with the histological structure of the corneal tissue, as observed by H&E staining (Fig. 1f). The expression of classic cell markers further validated the anatomical annotations of these sections (Fig. 1g and Additional file 2: Table S2). Using the matched scRNA-seq data as a reference, we mapped the 19 major cell clusters to our spatial transcriptomic data via the RCTD algorithm, assigning cell annotations to each spatial unit (bin 50) on the slides (Fig. 1h). The consequent mapping demonstrated notable concordance with identified anatomical regions and current biological knowledge, revealing spatial stratification and age-related cellular alterations in different corneal cell clusters (Fig. 1h).

Spatiotemporal changes in human corneal LSCs

LSCs include quiescent limbal stem cells (QLSCs) and active limbal stem cells (ALSCs). They are specialized cells located at the limbus of the cornea and are responsible for maintenance and repair following corneal injury. Cluster C9, characterized by its proliferative properties, represents transient amplifying cells (TACs) in the mitotic transition (Additional file 1: Fig. S2a). Cluster C10, characterized by increased expression of KRT15 and stem cell markers, were designated as LSCs. Cluster C11 was found to be predominantly located in the conjunctival epithelium in our spatial transcriptomic data and may represent the intermediate state of limbal cells differentiating into conjunctival epithelial cells (Additional file 1: Fig. S2b and c). Subsequent clustering of 2562 LSCs identified QLSCs and ALSCs, each with a distinct transcriptomic profile (Fig. 2a). Classic markers such as IFITM3 for QLSCs and ATF3 for ALSCs showed an exclusive expression pattern (Fig. 2b and c). GSEA revealed that the upregulated gene sets in the ALSCs were enriched in pathways related to “cell cycle”, “epithelial cell differentiation” and “muscle tissue development”, whereas the genes in the QLSCs were enriched in pathways related to “leukocyte chemotaxis”, “response to external stimulus” and “wound healing” (Fig. 2d).

Fig. 2
figure 2

Spatiotemporal changes in limbal stem cells (LSCs) during corneal aging. a UMAP visualization of two LSC subsets, active limbal stem cells (ALSCs) and quiescent limbal stem cells (QLSCs). b UMAP visualization showing the expression of ATF3 and IFITM3 in LSCs. c Violin plots showing the expression of ATF3 and IFITM3 in two LSC subsets. d GSEA enrichment plots showing the upregulated GO terms of ALSCs (top) and QLSCs (bottom). e UMAP visualization of three ALSC subsets. f Dot plot illustrating the expression of marker genes for each LSC subset. The color indicates the normalized expression level, while the dot size indicates the fraction of expressing cells. g Bar plot showing the results of the GO enrichment analysis of the marker genes of each LSC subset. h Line chart showing the proportion of different LSC subsets in each age stage. i Violin plot showing the stem cell differentiation score for each LSC subset. j Stacked plot showing the proportion of cell phases in each LSC subset. k Trajectory analysis of LSC subsets colored by pseudotime and subsets. l Spatial visualization of each LSC subset in each slide. m and n Immunofluorescence staining of KRT15 (red), CLDN4 (light blue), UBE2S (yellow), and SLC6A6 (green); nuclear DAPI staining is shown in blue

Full size image

We then focused on the ALSC subsets to further determine their unique roles in the intricate mechanisms governing their activation and differentiation. Further reclustering of 1068 ALSCs revealed three distinct subsets, namely, UBE2S+ ALSCs, CLDN4+ ALSCs and SLC6A6+ ALSCs, each with unique transcriptomic features (Fig. 2e and f). We projected these ALSC subsets back into the UMAP of all LSCs and visualized the expression patterns of their markers. Our analysis revealed that these markers exhibited distinct and specific expression profiles (Additional file 1: Fig. S3a–b). The SLC6A6+ ALSC subset, characterized by the expression of known limbal stem/progenitor cell markers such as SLC6A6, GPHA2, and KRT15, was associated with cellular respiration and mitochondrial function (Fig. 2f, g and Additional file 2: Table S3). The CLDN4+ ALSC subset, characterized by the expression of corneal epithelial markers (CLDN4, CCND1 and ELF3), was involved in processes related to intermediate filament organization and epithelial cell differentiation (Fig. 2f, g and Additional file 2: Table S3). The UBE2S+ ALSC subset, characterized by high expression of SOCS3, EGR1, and UBE2S, was enriched in genes related to cell proliferation and the response to epidermal growth factors (Fig. 2f, g and Additional file 2: Table S3). We observed intriguing changes in different LSC subsets during aging. In particular, we found that the QLSC had the highest proportion at 33Y and gradually decreased during aging. The populations of CLDN4+ ALSC and UBE2S+ ALSC significantly increased from 50Y and peaked at 68Y, whereas SLC6A6+ ALSC had the highest proportion in the advanced stages during aging (Fig. 2h).

Analysis of the stem cell differentiation and proliferation potential of each LSC subset revealed that UBE2S+ ALSCs exhibited enhanced differentiation capability and the highest proportion of actively proliferating cells, suggesting a more active nature (Fig. 2i and j). Trajectory analysis revealed a developmental trajectory showing that SLC6A6+ ALSCs and CLDN4+ ALSCs originated from QLSCs and evolved into UBE2S+ ALSCs (Fig. 2k), highlighting the dynamic nature of LSC state transitions. Spatial transcriptomic analysis and immunofluorescence staining confirmed the existence of these LSC subsets in the limbus (Fig. 2l–n). According to the immunofluorescence staining results, CLDN4+ ALSCs are located closer to the squamous epithelial layer, UBE2S+ ALSCs are located within the corneal epithelial layers, and SLC6A6+ ALSCs are distributed throughout the whole epithelial layer (Fig. 2m). Notably, the fluorescence intensity of SLC6A6+ ALSCs increased with increasing age (Fig. 2n).

We then used fuzzy clustering to identify transcriptomic changes associated with LSC aging and their functional characteristics. The optimal number of clusters was determined using Hubert's statistic value, which led to the identification of gene sets positively and negatively correlated with LSC aging (Additional file 1: Fig. S4a–b). Among these, the positively correlated gene sets were significantly enriched in functions such as “Translation”, “Positive regulation of gene expression” and “Regulation of cell cycle”, which may be related to the expansion of ALSCs in the aged groups (68Y and 88Y, Additional file 1: Fig. S4c). Conversely, the negatively correlated gene sets were associated with processes such as “Cell migration”, “Signal transduction” and immune-related pathways, including “Positive regulation of canonical NF-kappaB signal transduction” (Additional file 1: Fig. S4c). This pattern is consistent with the higher proportion of QLSCs observed in the younger groups.

We further investigated the molecular changes in the aging-related LSC subsets. Notably, “Pyruvate catabolic process” and “Pyruvate transmembrane transport” were negatively associated with QLSC aging (cluster 5, Additional file 1: Fig. S4 d–f). Previous studies have shown that the maintenance of quiescent stem cell homeostasis relies on mitochondrial pyruvate import [31]. These findings suggest that in younger samples, QLSCs may upregulate pyruvate metabolism processes to better maintain their functional state. During the aging process of SLC6A6-ALSCs, intrinsic functions such as “Mitochondrial translation” were further activated, while apoptotic processes were also enriched, which may correlate with the aging of SLC6A6-ALSCs (Additional file 1: Fig. S4 g–i).

We further investigated the key TFs that regulate the ultimate transition of LSCs via the pySCENIC algorithm (Fig. 3a). Our analysis revealed several TFs that were uniquely enriched in different subsets of LSCs. In particular, several TFs known to promote differentiation are expressed predominantly in UBE2S+ ALSC. In addition, the expression of hypoxia inducible factor 1 alpha (HIF1A), a TF critical for wound healing, was significantly increased in QLSCs (Fig. 3a). We also discovered a HIF1A-associated regulatory module involved in the regulation of responses to interferon-beta, hypoxia, angiogenesis, and wound healing (Fig. 3b). The localized pattern of HIF1A expression in QLSCs led to the subdivision of these cells into HIF1A-positive (HIF1A+ QLSCs) and HIF1A-negative (HIF1A− QLSCs) subsets (Fig. 3c). There was little difference in the expression of the QLSC marker IFITM3 between these two subsets (Fig. 3d). However, HIF1A+ QLSCs showed significantly enhanced wound healing and leukocyte chemotaxis abilities (Fig. 3d).

Fig. 3
figure 3

Spatiotemporal changes in quiescent limbal stem cells (QLSCs) during corneal aging. a Heatmap showing the expression of transcription factors (TFs) for each LSC subset. b Heatmap showing the results of the GO enrichment analysis of HIF1A and its downstream target genes. c UMAP visualization showing the expression of HIF1A in LSCs. d Violin plots showing the expression of IFITM3 and the signature score of wound healing and leukocyte chemotaxis in HIF1A− QLSCs and HIF1A+ QLSCs. e GSEA enrichment plot showing the upregulated GO terms of HIF1A+ QLSCs. f Chord diagram showing the interaction between HIF1A+ QLSCs, HIF1A− QLSCs and vascular endothelial cells. g Lollipop plot showing differences in interaction strength between each cell type and between HIF1A+ QLSCs and HIF1A− QLSCs. h Spatial visualization of HIF1A− QLSCs and HIF1A+ QLSCs, with insets showing enlarged areas of the selected regions. i Immunofluorescence staining of KRT15 (green), HIF1A (red), and CCL14 (lake blue), with nuclei stained with DAPI (blue)

Full size image

GSEA of HIF1A+ QLSCs and HIF1A− QLSCs revealed “positive regulation of blood vessel endothelial cell migration” as the most significantly enriched GO term of the HIF1A+ QLSCs (Fig. 3e). By inferring cell–cell communication, we showed that compared with HIF1A− QLSCs, HIF1A+ QLSCs exhibited increased crosstalk with vascular endothelial cells (Fig. 3f), with vascular endothelial cells showing the most pronounced difference in communication strength between the two cell types (Fig. 3g). To contextualize these findings within the tissue architecture, we analyzed the spatial localization of HIF1A+ QLSCs, HIF1A− QLSCs and vascular endothelial cells within the cornea and found that HIF1A+ QLSCs were predominantly located close to the limbal epithelium and spatially adjacent to vascular endothelial cells, whereas HIF1A− QLSCs were more dispersed in our spatial transcriptomics data (Fig. 3h). Immunofluorescence staining further confirmed the proximity of HIF1A+ QLSCs to the corneal epithelium and near vascular endothelial cells (CCL14), revealing their specialized niche within the limbal area (Fig. 3i). Based on the enhanced interaction between HIF1A+ QLSCs and vascular endothelial cells and the specific localization of HIF1A+ QLSCs within the limbal niche, our results delineate the spatial distribution of HIF1A+ QLSCs in relation to the limbal epithelium and detail their intricate interactions with neighboring vascular endothelial cells.

Diversity and spatial location of corneal basal cells

To elucidate the spatiotemporal dynamics of corneal epithelial cells, we reanalyzed five identified subsets of corneal epithelial cells in a UMAP plot (Fig. 4a). Clusters C6, C7 and C8, characterized by high expression of GJA1, COL17A1, and GJB2, which are markers typically found in the basal epithelial cells, were classified into distinct corneal basal epithelial subsets (Fig. 4bc, Additional file 1: Fig. S5a). KRT3, a marker of the suprabasal corneal epithelium, was predominantly expressed in C4, indicating that this subset comprises wing cells (also known as corneal suprabasal epithelial cells) (Fig. 4b and c). Cluster C5, identified by the presence of KRT24, was classified as corneal squamous epithelial cells (Fig. 4b and c). Our analysis revealed the developmental progression of corneal epithelial cells, showing their differentiation from basal epithelial cells to wing cells, culminating in their maturation into squamous epithelial cells (Fig. 4d).

Fig. 4
figure 4

Diversity and spatial location of corneal basal cells during aging. a UMAP visualization of five corneal epithelial cell (CEC) clusters. b UMAP visualization showing the expression of GJA1, KRT24 and KRT3 in CECs. c Violin plots showing the expression of GJA1, KRT24 and KRT3 in each CEC subset. d Trajectory analysis of CEC subsets colored by pseudotime and subsets. e Dot plot illustrating the expression of marker genes for three basal subsets (left). The color indicates the normalized expression level, while the dot size indicates the fraction of expressing cells. Bar plot displaying the results of the GO enrichment analysis of the marker genes (right). f Spatial visualization of each basal subset in each slide, with insets showing enlarged areas of the selected regions. g Immunofluorescence staining of NTRK2 (purple) and DAPI (blue) in corneal sections. h Schematic representation of the corneal regions. i Western blot showing NTRK2 protein levels in young (8 weeks old) and aged (14 months old) mouse corneas. j Histogram showing the comparative expression levels of NTRK2 in the corneas of young and aged mice. The error bars indicate the standard error of the mean (p* < 0.05). k Violin plot showing the epithelial cell differentiation score for each basal subset. Chord diagram showing the interaction between three basal subsets, wing cells, and squamous cells. l GSEA enrichment plot showing the enriched GO terms associated with genes whose expression was upregulated in NTRK2+THBS1+ basal cells compared to that in NTRK2+ THBS1− basal cells. m Spatial visualization of NTRK2+THBS1 + basal cells, NTRK2+THBS1− basal cells, and other CECs in each slide, with insets showing enlarged areas of the selected regions

Full size image

We further investigated the transcriptomic heterogeneity and functional diversity of these corneal basal subsets (Additional file 2: Table S4). Specifically, basal C6 cells expressed high levels of the stem cell marker KRT15, together with enrichment in pathways associated with epithelial cell differentiation and cell differentiation (Fig. 4e, Additional file 1: Fig. S5b). Basal C8 cells are characterized by increased expression of lymphocyte antigen 6 family member D (LY6D) and the inflammatory mediators S100A8 and S100A9, which are involved in processes such as “neutrophil aggregation”, the “innate immune response”, and “leukocyte migration involved in the inflammatory response” (Fig. 4e, Additional file 1: Fig. S5b). Basal C7 cells, marked by the putative limbal progenitor cell marker S100A2, chemokine ligand 14 (CXCL14) and the corneal basal epithelial marker NTRK2, were associated with pathways related to “aging”, “cell adhesion”, and “apoptotic process” (Fig. 4e, Additional file 1: Fig. S5b). These three basal subsets were classified as KRT15+ (C6), NTRK2+ (C7), and LY6D+ (C8) basal cells based on their unique molecular signature (Fig. 4e).

We observed an overall increase in the proportion of squamous cells with aging (Additional file 1: Fig. S6a). By examining the molecular changes during aging, we found that squamous epithelial cells showed a significant enrichment in ATP production and energy metabolism functions with aging, including “Proton motive force-driven mitochondrial ATP synthesis”, “Aerobic respiration”, and “Response to oxidative stress”. This suggests that the gradual increase in corneal squamous epithelial cells during aging is associated with the activation of energy metabolism pathways (cluster 5, Additional file 1: Fig. S6b–d). Notably, the three basal subsets showed distinct enrichment patterns during the aging process. LY6D+ basal cells were predominantly enriched at 33Y, KRT15+ basal cells became dominant at 50 and 68Y, whereas the proportion of the NTRK2+ basal cells increased sharply at 88Y (Additional file 1: Fig. S6a).

The spatial distribution of the basal cell subsets was clearly stratified, with a prevalence of basal cells in the peripheral cornea at 33 years of age, which transitioned to central accumulation in older samples in the spatial transcriptomics data (Fig. 4f and Additional file 1: Fig. S7a, b). Specifically, NTRK2+ basal cells appeared to dominate the central cornea in aging corneas, whereas KRT15+ basal cells were predominantly found in the periphery (Fig. 4f and Additional file 1: Fig. S7c). Consistent with the spatial transcriptomics findings, immunofluorescence staining revealed that NTRK2+ basal cells persisted in the peripheral corneal epithelium across all ages, but an increasing trend with age was observed in the central corneal epithelium (Fig. 4g). Building on our findings in human corneas, we similarly assessed NTRK2 expression across different corneal regions in young and aged mouse corneas (Fig. 4h). Western blot results revealed increased NTRK2 protein levels in aged mice, which were marked by a transition from peripheral prominence in young mice to increased central expression in aged mice (Fig. 4i and j). The histogram shows a statistically significant increase in NTRK2 expression with age (Fig. 4j), suggesting that the changes in NTRK2 expression during corneal aging may follow a similar pathway in mice as that observed in humans. GO enrichment analysis of the DEGs revealed distinct functional roles of basal cells in different corneal regions in our spatial transcriptomics data, with central corneal basal cells primarily involved in defense and response to viral infection and peripheral corneal basal cells involved in intermediate filament organization and keratinization (Additional file 1: Fig. S7 d).

To further elucidate the multicellular communication networks and mutual regulatory interactions within the corneal epithelium, we applied a classical approach to calculate CNs on shared tissue niches in our spatial transcriptomics data. According to the canonical approach, we defined the window of each individual spot as the N = 10 nearest spatial neighbors. Next, the optimal number of CN clusters k for each sample was determined based on the highest average silhouette score (Additional file 1: Fig. S8a). We identified 8 CNs in the 33Y sample, 5 CNs in the 50Y sample, 7 CNs in the 68Y sample, and 6 CNs in the 88Y sample (Additional file 1: Fig. S8b). Through correlation analysis and anatomical region enrichment analysis, we identified several major CNs associated with different corneal regions, including corneal epithelium-related CNs (33Y-CN1, 68Y-CN6, 50Y-CN3, 88Y-CN2), corneal endothelium-related CNs (33Y-CN8, 68Y-CN7, 88Y-CN4), and corneal stroma-related CNs (33Y-CN6, 50Y-CN4, 68Y-CN2, 88Y-CN1) (Additional file 1: Fig. S8c). In corneal epithelium-related CNs (Additional file 1: Fig. S8 d), we observed that the 33Y-CN1 predominantly consisted of wing cells and corneal squamous epithelial cells, whereas 50Y-CN3, 68Y-CN6, and 88Y-CN2 included corneal basal cells, wing cells, and corneal squamous epithelial cells. This is consistent with our observations that in the 33Y sample, corneal basal cells were mainly located in the peripheral cornea, further away from the mature corneal epithelium (wing cells and corneal squamous epithelial cells). Cell–cell communication analysis further revealed that the interaction strength between corneal basal cells and wing cells/corneal squamous epithelial cells was the lowest in the 33Y sample (Additional file 1: Fig. S8e). In addition, the corneal epithelium-related CNs in the 88Y sample, 88Y-CN2, primarily consists of all subtypes of corneal epithelial cells and also SLC6A6-ALSCs, which is consistent with our findings in the single-cell analysis, showing that SLC6A6-ALSCs was significantly upregulated in the elderly samples. These results suggest that in young corneal samples, mature corneal epithelial cells (wing cells and corneal squamous epithelial cells) and immature corneal epithelial cells (corneal basal cells) form independent spatial CNs. In aging corneas, the upregulation of SLC6A6-ALSCs and its close neighborhood relationship with other corneal epithelial cells may play a key regulatory role of SLC6A6-ALSCs during the aging process.

Furthermore, the epithelial differentiation score of NTRK2+ basal cells was lower than that of KRT15+ and LY6D+ basal cells, and NTRK2+ basal cells exhibited weaker communication with differentiated corneal epithelial cells (wing cells and squamous cells) (Fig. 4k). Additional analysis of receptor‒ligand interactions between basal cells and differentiated corneal epithelial cells revealed a unique cell‒cell interaction pattern in NTRK2 + basal cells (Additional file 1: Fig. S9a). Notably, THBS1-mediated interactions had the strongest effect on NTRK2+ basal cells in 88-year-old corneas (Additional file 1: Fig. S9b). THBS1 was predominantly expressed in NTRK2 + basal cells and served as a hub gene in the THBS1-associated gene module, which was correlated with cell adhesion and migration (Additional file 1: Fig. S9c–e). Differential gene expression analysis of NTRK2+THBS1+ basal cells and NTRK2+THBS1− basal cells revealed significant upregulation of genes associated with “epithelial cell differentiation” and “regulation of cell adhesion” in NTRK2+ THBS1+ basal cells (Fig. 4l). These results suggest that THBS1 expression in NTRK2 + basal cells may be crucial for maintaining their functionality and differentiation capacity through specific cell‒cell interactions. Further investigation into the influence of THBS1 on the spatial distribution of NTRK2+ basal cells via spatial transcriptomics data revealed a layered distribution pattern between NTRK2+THBS1− and NTRK2+ THBS1+ basal cells. The latter were closer to differentiated corneal epithelial cells in all the samples (Fig. 4m), suggesting that THBS1 may be involved in mediating the adhesive, migratory and differentiation functions of NTRK2+ basal cells.

Characterization of stromal cell subsets

We reclustered the stromal cells into five distinct subsets, excluding a mixed subset of Schwann cells identified within this population (Fig. 5a and b). All identified stromal subsets were characterized by high expression of the stromal cell marker decorin (DCN) (Fig. 5b). DEG analysis of these stromal subsets (Additional file 2: Table S5) revealed that a subset of C4 cells were corneal stromal stem cells (CSSCs), characterized by upregulation of the stem cell markers KRT15, S100A2, and PAX6 (Fig. 5c). CSSCs were predominantly distributed along the corneal periphery, and their number was significantly higher in young corneas than in aged samples (68 years and 88 years) (Fig. 5e). The remaining three stromal cell subsets showed unique characteristics and functional profiles. Stromal-1, characterized by the upregulation of the fibroblast marker CFD and the collagen-encoded genes COL14A1 and COL1A2, was involved in “collagen fibril organization” and “extracellular matrix organization” (Fig. 5c and d). Spatial transcriptomics data indicated that these cells were primarily located within the scleral stroma; therefore, they were classified as scleral fibroblasts (Fig. 5e). Stromal-2 cells, distinguished by the expression of SAA1 and KERA, were identified as corneal stromal cells and were distributed throughout the corneal stroma in the spatial transcriptomic data (Fig. 5c and e). GO enrichment analysis revealed that Stromal- 2 was associated with “cell adhesion” and “regulation of tissue remodeling,” which is consistent with the known characteristics of corneal stromal cells (Fig. 5d). Stromal-3, which expresses RBP4, NR2F1, MYOC, and DDIT4, was predominantly localized to the ciliary body (Fig. 5e) and was associated with “negative regulation of cell proliferation”, “negative regulation of signal transduction”, and “anatomical structure development” (Fig. 5c and d). Leveraging spatial information, we clearly delineated scleral stromal cells, corneal stromal cells, and ciliary body stromal cells, and identified their respective specific markers as well as their unique features and functional profiles.

Fig. 5
figure 5

Characterization of stromal cells during aging. a UMAP visualization of five stromal cell subsets. b UMAP visualization showing the expression of DCN and S100B in stromal cells. c Dot plot illustrating the expression of marker genes for four stromal cell subsets and corneal stromal stem cells (CSSCs). The color indicates the normalized expression level, while the dot size indicates the proportion of expressing cells. d Heatmap showing the expression of differentially expressed genes (DEGs) for three stromal cell subsets (left). Bar plot showing the results of GO enrichment analysis of the DEGs of three stromal cell subsets (right). e Spatial visualization of each stromal cell subset in each slide. f Heatmap showing the expression of collagen, extracellular matrix (ECM) proteoglycans and elastic fiber-encoded genes for the four stromal cell subsets in the spatial transcriptomics data. g Heatmap showing the expression of collagen, ECM proteoglycans and elastic fiber-encoded genes for stromal-1 (left) and stromal-2 (right) across all age groups in the spatial transcriptomics data

Full size image

Given the critical role of stromal cells in the generation of a collagen-rich ECM, we further investigated the variation in essential stromal functions between different subsets of stromal cells and their changes with age. The major collagen types, collagen type I (COL1A1 and COL1A2) and type V (COL5A1 and COL5A2), were predominantly expressed by scleral stromal cells and corneal stromal cells, respectively (Fig. 5f). ECM proteoglycans, which are critical for the structural integrity and function of the stromal ECM, were significantly upregulated in corneal stromal cells (Fig. 5f). Furthermore, spatial transcriptomics data revealed a significant increase in the expression of collagen type I, collagen type V, and ECM in both scleral stromal cells and corneal stromal cells from aged samples (68 and 88 years) (Fig. 5g).

Spatial, molecular and functional heterogeneity of corneal endothelial cells

To investigate potential spatial-functional heterogeneity, we categorized the corneal endothelium into peripheral and central regions based on its spatial distribution in our spatial transcriptomics data (Fig. 6a). The corneal endothelium showed the highest level of mitochondrial energy metabolism in the cornea (Fig. 6b), with the peripheral endothelium showing a higher proportion of mitochondrial gene expression (Fig. 6c). Differential gene expression analysis of the peripheral and central endothelium revealed genes enriched in “cellular respiration”, “ATP synthesis coupled proton transport”, and “response to calcium ions,” whereas the central endothelium was involved in processes such as “negative regulation of glycolytic process” (Fig. 6d and e). This finding suggests that there are differences in metabolic capacity between the peripheral and central endothelium. We further collected well-defined metabolic pathways from KEGG to investigate the differences in metabolic patterns between the two classes of endothelium. Only oxidative phosphorylation was significantly upregulated in the peripheral endothelium (Fig. 6f), confirming that oxidative phosphorylation is the major function influencing their metabolic differences.

Fig. 6
figure 6

Spatial, molecular and functional heterogeneity of corneal endothelial cells. a Spatial visualization delineating the peripheral and central corneal endothelium. b Spatial visualization of endothelial cells with different percentages of mitochondrial gene expression. c Violin plots showing the percentages of mitochondrial gene expression in the peripheral and central corneal endothelium. d Scatter plot showing log2-fold changes in gene expression in the peripheral and corneal endothelium. e Bar plot showing the results of the GO enrichment analysis of the DEGs associated with the peripheral and central corneal endothelium. f GSEA enrichment plot displaying the significantly upregulated metabolic pathways of the peripheral corneal endothelium. g Lollipop chart showing the functional scores of the peripheral and central corneal endothelium across different age groups. h Corneal endothelial cells costained with CA3 (red) and the Na–K-ATPase subunit (green), with DAPI in blue (left). A statistical image of the Na–K-ATPase immunofluorescence results is shown on the right (p*** < 0.001). i Immunofluorescence staining of mouse corneal whole mounts. NDUFV2 (I) (red) and SDHB (II) (lake blue) are expressed at higher levels in the peripheral corneal endothelium than in the central region. An immunofluorescence statistical diagram of NDUFV2 (I) and SDHB (II) is shown on the right (p*** < 0.001)

Full size image

Further functional scoring revealed that the peripheral endothelium exhibited pronounced ion transport activity, oxidative phosphorylation, and ATP energy metabolism compared to the central region (Fig. 6g). This difference suggests that the peripheral endothelium is more metabolically active than the central endothelium. However, the functional differences between these two regions decreased with increasing age (Fig. 6g). Immunostaining for Na–K-ATPase, which is used to assess ion pump performance, corroborated these spatial transcriptomic findings by showing higher expression in the corneal endothelium on both sides than in the central endothelium (Fig. 6h). Immunofluorescence staining of 8-week-old mouse corneal whole mounts revealed higher expression levels of NDUFV2 and SDHB in the peripheral corneal endothelium than in the central region (Fig. 6i), indicating a greater capacity for oxidative phosphorylation. This finding is consistent with the spatial transcriptomics data showing increased metabolic activity in the peripheral corneal endothelium. These results highlight the functional characteristics of different regions of the corneal endothelium.

By analyzing transcriptional changes in the corneal endothelium with aging, we identified two gene sets negatively correlated with corneal endothelial aging (Additional file 1: Fig. S10a–b). Functional enrichment analysis of the four selected gene sets revealed that endothelial cell aging is associated with processes such as metabolism and the cell cycle, whereas processes such as “Translation”, “Intracellular zinc ion homeostasis”, and “Response to manganese ion” were found to be gradually downregulated during endothelial cell aging (Additional file 1: Fig. S10c).

Discussion

Although bulk and single-cell omics studies in animal models, corneal organoids or human corneal limbal tissue have provided some insight into corneal biology [7, 32,33,34], the understanding of the spatial, cellular and molecular architecture of the human cornea remains incomplete due to the limited availability of intact human corneas across a broad age spectrum. This study addresses this critical gap by providing a comprehensive, high-resolution spatial transcriptomic and cellular atlas of the adult whole human cornea. By incorporating scRNA-seq and scStereo-seq data, this spatiotemporal atlas reveals the spatial distribution of major cell types and their subclusters, as well as specific gene expression signatures of the healthy human cornea, thereby providing a valuable resource for advancing research in corneal pathology, transplantation, senescence and regenerative medicine.

Regeneration and homeostasis of the corneal epithelium depend on the presence of LSCs [35, 36], whereas limbal stem cell deficiency can ultimately lead to corneal blindness [37]. Notably, the corneal limbus has been found to be organized into two distinct stem cell compartments [38]. Previous studies have identified two distinct populations of stem/progenitor cells in the corneal limbus: QLSCs at the peripheral corneal limbus, which play a role in injury repair, and ALSCs in the inner corneal limbus, which are involved in corneal epithelial regeneration [7, 39]. Given the critical role of QLSCs as reservoirs for tissue regeneration and their significant contribution to boundary maintenance, understanding their localization and regulation is of paramount importance for advancing the field of corneal biology and developing therapeutic strategies for the treatment of corneal diseases. Our results confirm the existence and major function of these two LSC populations in the human cornea and reveal that aging significantly affects the function and composition of these stem cell populations, disrupting tissue homeostasis and compromising regenerative capacity [40, 41]. Previous studies have shown that the number and activation of quiescent stem cells are disrupted during aging, leading to impaired regeneration in multiple tissues [40, 42,43,44,45]. Consistent with these findings, our study revealed a significant reduction in the number of QLSCs in aged samples. Based on our further observations of QLSCs, we identified HIF1A as a key TF in QLSCs associated with enhanced wound healing and leukocyte chemotaxis. Importantly, our research revealed that a specialized subset of HIF1A-positive QLSCs, strategically located at the limbal periphery near vascular structures strongly interacts with vascular endothelial cells and has an enhanced ability to regulate endothelial migration. Deficiency or dysfunction of LSCs, as observed in conditions such as aniridia, Stevens-Johnson syndrome or corneal chemical burn, has been identified as a significant factor contributing to corneal neovascularization [46, 47]. We suggest that these QLSCs are critical for maintaining the avascular integrity of the cornea. These findings may pave the way for new therapeutic strategies to help maintain corneal transparency and prevent neovascularization, both of which are critical for maintaining visual acuity and preventing blindness. A better understanding of the molecular biology of LSC regulation, particularly in relation to age-related changes in LSCs, and improvements in LSC bioengineering will be paramount in improving the clinical outcomes of LSC treatments.

Our results also challenge the prevailing assumptions about the uniformity of the corneal epithelial layer distribution with age, revealing a more complex, age-dependent migration pattern of basal cells from the periphery to the central cornea [48,49,50]. By defining three distinct basal cell subsets-KRT15+ basal, NTRK2+ basal, and LY6D+ basal, we demonstrated that basal cells of the epithelium migrate from the periphery to the central cornea and that NTRK2+ basal cells become the predominant subset in the central cornea with increasing age. Recent studies have also shown that corneal epithelial cells typically originate from the peripheral region of the cornea and gradually migrate toward the center (“centripetal pattern”) with each epithelial renewal [7, 51, 52]. Based on these findings, we propose the “spatiotemporal centripetal pattern” as a new paradigm to explain the age-related centripetal migration of corneal epithelial cells, which supports the maintenance of corneal homeostasis through the dynamic transformation of the LSC subpopulation and epithelial cells. We also investigated the key role of THBS1 in NTRK2+ basal cells, which is related to cell migration and apoptosis. Interestingly, NTRK2+ THBS1+ basal cells were found in close proximity to differentiated corneal epithelial cells. Recent studies have elucidated the central role of THBS1 in the aging process. Porpiglia et al. reported that THBS1 interacts with CD47 on aging muscle stem cells, suppressing their proliferative capacity and impairing their muscle regeneration potential [53]. These findings highlight the importance of THBS1 in cell senescence and provide new directions for future therapeutic strategies. We hypothesize that the modulation of THBS1 expression in NTRK2+ basal cells may influence the regenerative potential of the corneal epithelium.

This study has several limitations. First, the sample size was relatively small and limited to male subjects. Future studies should be expanded to include larger and more diverse samples to strengthen these findings and to investigate potential sex differences in corneal aging. Second, although scRNA-seq provides transcriptomic information at single-cell resolution, the relatively limited sequencing depth did not capture the whole transcriptome, which may affect cell annotation and state interpretation. Therefore, for the critical cell populations identified in our study, we performed additional validation using immunofluorescence to confirm their presence in corneal tissue. Finally, our analyses were limited to corneal sections along the nasotemporal axis. Further investigation of different regions of the whole cornea is needed to gain a more comprehensive understanding of regional differences in corneal aging.

Conclusions

In conclusion, our comprehensive analysis of the single-cell and spatial transcriptomic landscape of the human cornea at different ages provides a high-resolution spatiotemporal atlas that reveals significant differences in gene expression, cellular composition and spatial organization of major corneal cell types during aging. This study provides the first reference atlas of aging hallmarks in the whole human cornea, advancing our understanding of the cellular and molecular dynamics of human corneal aging and shedding light on the molecular mechanisms underlying age-related eye diseases.

Data availability

The scRNA-seq and scStereo-seq data generated in this study have been deposited in the Genome Sequence Archive for Human (GSA-Human) repository under the accession number HRA007528 (https://ngdc.cncb.ac.cn/gsa-human/browse/HRA007528) [54] and HRA007620 (https://ngdc.cncb.ac.cn/gsa-human/browse/HRA007620) [55].

Abbreviations

ALSCs:

Active limbal stem cells

ECM:

Extracellular matrix

GJB4:

Gap junction protein beta 4

GSEA:

Gene set enrichment analysis

H&E:

Hematoxylin and eosin

HIF1A:

Hypoxia-inducible factor 1-alpha subunit

HES1:

Hes family bHLH transcription factor 1

IFITM3:

Interferon-induced transmembrane protein 3

JAG1:

Jagged 1 (delta-like 1 homolog)

KERA:

Keratocan

KRT15:

Keratin 15

KRT3:

Keratin 3

KRT24:

Keratin 24

LY6D:

Lymphocyte antigen 6 complex, locus D

LSCs:

Limbal stem cells

LSCD:

Limbal stem cell deficiency

MSigDB:

Molecular Signatures Database

NTRK2:

Neurotrophic tyrosine kinase receptor type 2

NDUFV2:

NADH dehydrogenase [ubiquinone] flavoprotein 2, mitochondrial

NR2 F1:

Nuclear receptor subfamily 2, group F, member 1

PAX6:

Paired box 6

PPI:

Protein-protein interaction

QLSC:

Quiescent limbal stem cells

RBP4:

Retinol binding protein 4, plasma

SDHB:

Succinate dehydrogenase complex, subunit B, integral membrane protein

SAA1:

Serum amyloid A1

SCENIC:

Single-Cell ENrichment of Gene Set and Inference of Gene Regulatory Network

scStereo-seq:

Single-cell SpaTial Enhanced REsolution Omics-sequencing

scRNA-seq:

Single-cell RNA sequencing

SLC6A6:

Solute carrier family 6, member 6

SOCS3:

Suppressor of cytokine signaling 3

S100A2:

S100 calcium binding protein A2

S100A8:

S100 calcium binding protein A8

S100A9:

S100 calcium binding protein A9

TACs:

Transient amplifying cells

THBS1:

Thrombospondin 1

UBE2S:

Ubiquitin conjugating enzyme E2S

References

  1. Li S, Zenkel M, Kruse FE, Gießl A, Schlötzer-Schrehardt U. Identification, isolation, and characterization of melanocyte precursor cells in the human limbal stroma. Int J Mol Sci. 2022;23(7): 3756.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  2. Sheerin AN, Smith SK, Jennert-Burston K, Brook AJ, Allen MC, Ibrahim B, et al. Characterization of cellular senescence mechanisms in human corneal endothelial cells. Aging Cell. 2012;11(2):234–40.

    Article  CAS  PubMed  Google Scholar 

  3. Sagga N, Kuffová L, Vargesson N, Erskine L, Collinson JM. Limbal epithelial stem cell activity and corneal epithelial cell cycle parameters in adult and aging mice. Stem cell research. 2018;33:185–98.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Hillenaar T, van Cleynenbreugel H, Remeijer L. How normal is the transparent cornea? Effects of aging on corneal morphology. Ophthalmology. 2012;119(2):241–8.

    Article  PubMed  Google Scholar 

  5. Ljubimov AV, Saghizadeh M. Progress in corneal wound healing. Prog Retin Eye Res. 2015;49:17–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. Price MO, Mehta JS, Jurkunas UV, Price FW Jr. Corneal endothelial dysfunction: Evolving understanding and treatment options. Prog Retin Eye Res. 2021;82:100904.

    Article  PubMed  Google Scholar 

  7. Altshuler A, Amitai-Lange A, Tarazi N, Dey S, Strinkovsky L, Hadad-Porat S, et al. Discrete limbal epithelial stem cell populations mediate corneal homeostasis and wound healing. Cell Stem Cell. 2021;28(7):1248-61.e8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Lin JB, Shen X, Pfeifer CW, Shiau F, Santeford A, Ruzycki PA, et al. Dry eye disease in mice activates adaptive corneal epithelial regeneration distinct from constitutive renewal in homeostasis. Proc Natl Acad Sci U S A. 2023;120(2): e2204134120.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  9. Li M, Guo H, Wang B, Han Z, Wu S, Liu J, et al. The single-cell transcriptomic atlas and RORA-mediated 3D epigenomic remodeling in driving corneal epithelial differentiation. Nat Commun. 2024;15(1):256.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Lein E, Borm LE, Linnarsson S. The promise of spatial transcriptomics for neuroscience in the era of molecular cell typing. Science. 2017;358(6359):64–9.

    Article  CAS  PubMed  Google Scholar 

  11. Crosetto N, Bienko M, van Oudenaarden A. Spatially resolved transcriptomics and beyond. Nat Rev Genet. 2015;16(1):57–66.

    Article  CAS  PubMed  Google Scholar 

  12. Chen A, Sun Y, Lei Y, Li C, Liao S, Meng J, et al. Single-cell spatial transcriptome reveals cell-type organization in the macaque cortex. Cell. 2023;186(17):3726-43 e24.

    Article  CAS  PubMed  Google Scholar 

  13. Chen A, Liao S, Cheng M, Ma K, Wu L, Lai Y, et al. Spatiotemporal transcriptomic atlas of mouse organogenesis using DNA nanoball-patterned arrays. Cell. 2022;185(10):1777-92.e21.

    Article  CAS  PubMed  Google Scholar 

  14. Wei X, Fu S, Li H, Liu Y, Wang S, Feng W, et al. Single-cell Stereo-seq reveals induced progenitor cells involved in axolotl brain regeneration. Science. 2022;377(6610):eabp9444.

    Article  CAS  PubMed  Google Scholar 

  15. Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573-87.e29.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Wolock SL, Lopez R, Klein AM. Scrublet: computational identification of cell doublets in single-cell transcriptomic data. Cell Syst. 2019;8(4):281-91.e9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Korsunsky I, Millard N, Fan J, Slowikowski K, Zhang F, Wei K, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16(12):1289–96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Sherman BT, Hao M, Qiu J, Jiao X, Baseler MW, Lane HC, et al. DAVID: a web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res. 2022;50(W1):W216–21.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. Kumar L, Futschik ME. Mfuzz: a software package for soft clustering of microarray data. Bioinformation. 2007;2(1):5–7.

    Article  PubMed  PubMed Central  Google Scholar 

  20. Schurch CM, Bhate SS, Barlow GL, Phillips DJ, Noti L, Zlobec I, et al. Coordinated cellular neighborhoods orchestrate antitumoral immunity at the colorectal cancer invasive front. Cell. 2020;182(5):1341-59 e19.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Qiu X, Mao Q, Tang Y, Wang L, Chawla R, Pliner HA, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. 2017;14(10):979–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  22. Van de Sande B, Flerin C, Davie K, De Waegeneer M, Hulselmans G, Aibar S, et al. A scalable SCENIC workflow for single-cell gene regulatory network analysis. Nat Protoc. 2020;15(7):2247–76.

    Article  PubMed  Google Scholar 

  23. Kumar N, Mishra B, Athar M, Mukhtar S. Inference of gene regulatory network from single-cell transcriptomic data using pySCENIC. Methods Mol Biol. 2021;2328:171–82.

    Article  CAS  PubMed  Google Scholar 

  24. Aibar S, Gonzalez-Blas CB, Moerman T, Huynh-Thu VA, Imrichova H, Hulselmans G, et al. SCENIC: single-cell regulatory network inference and clustering. Nat Methods. 2017;14(11):1083–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Jin S, Guerrero-Juarez CF, Zhang L, Chang I, Ramos R, Kuan CH, et al. Inference and analysis of cell-cell communication using Cell Chat. Nat Commun. 2021;12(1):1088.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Gong C, Li S, Wang L, Zhao F, Fang S, Yuan D, et al. SAW: an efficient and accurate data analysis workflow for Stereo-seq spatial transcriptomics. GigaByte. 2024;2024:gigabyte111.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21.

    Article  CAS  PubMed  Google Scholar 

  28. Cable DM, Murray E, Zou LS, Goeva A, Macosko EZ, Chen F, et al. Robust decomposition of cell type mixtures in spatial transcriptomics. Nat Biotechnol. 2022;40(4):517–26.

    Article  CAS  PubMed  Google Scholar 

  29. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdottir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27(12):1739–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Maiti G, Monteiro de Barros MR, Hu N, Dolgalev I, Roshan M, Foster JW, et al. Single cell RNA-seq of human cornea organoids identifies cell fates of a developing immature cornea. PNAS Nexus. 2022;1(5):pgac246.

    Article  PubMed  PubMed Central  Google Scholar 

  31. Petrelli F, Scandella V, Montessuit S, Zamboni N, Martinou JC, Knobloch M. Mitochondrial pyruvate metabolism regulates the activation of quiescent adult neural stem cells. Sci Adv. 2023;9(9):eadd5220.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  32. Swarup A, Phansalkar R, Morri M, Agarwal A, Subramaniam V, Li B, et al. Single-cell transcriptomic analysis of corneal organoids during development. Stem Cell Reports. 2023;18(12):2482–97.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  33. Collin J, Queen R, Zerti D, Bojic S, Dorgau B, Moyse N, et al. A single cell atlas of human cornea that defines its development, limbal progenitor cells and their interactions with the immune cells. Ocul Surf. 2021;21:279–98.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Li DQ, Kim S, Li JM, Gao Q, Choi J, Bian F, et al. Single-cell transcriptomics identifies limbal stem cell population and cell types mapping its differentiation trajectory in limbal basal epithelium of human cornea. Ocul Surf. 2021;20:20–32.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Davanger M, Evensen A. Role of the pericorneal papillary structure in renewal of corneal epithelium. Nature. 1971;229(5286):560–1.

    Article  CAS  PubMed  Google Scholar 

  36. Cotsarelis G, Cheng SZ, Dong G, Sun TT, Lavker RM. Existence of slow-cycling limbal epithelial basal cells that can be preferentially stimulated to proliferate: implications on epithelial stem cells. Cell. 1989;57(2):201–9.

    Article  CAS  PubMed  Google Scholar 

  37. Majo F, Rochat A, Nicolas M, Jaoude GA, Barrandon Y. Oligopotent stem cells are distributed throughout the mammalian ocular surface. Nature. 2008;456(7219):250–4.

    Article  CAS  PubMed  Google Scholar 

  38. Farrelly O, Suzuki-Horiuchi Y, Brewster M, Kuri P, Huang S, Rice G, et al. Two-photon live imaging of single corneal stem cells reveals compartmentalized organization of the limbal niche. Cell Stem Cell. 2021;28(7):1233-47.e4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  39. Ouyang H, Xue Y, Lin Y, Zhang X, Xi L, Patel S, et al. WNT7A and PAX6 define corneal epithelium homeostasis and pathogenesis. Nature. 2014;511(7509):358–61.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  40. Goodell MA, Rando TA. Stem cells and healthy aging. Science. 2015;350(6265):1199–204.

    Article  CAS  PubMed  Google Scholar 

  41. Brunet A, Goodell MA, Rando TA. Ageing and rejuvenation of tissue stem cells and their niches. Nat Rev Mol Cell Biol. 2023;24(1):45–62.

    Article  CAS  PubMed  Google Scholar 

  42. de Morree A, Rando TA. Regulation of adult stem cell quiescence and its functions in the maintenance of tissue integrity. Nat Rev Mol Cell Biol. 2023;24(5):334–54.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Ji J, Ho BS, Qian G, Xie XM, Bigliardi PL, Bigliardi-Qi M. Aging in hair follicle stem cells and niche microenvironment. J Dermatol. 2017;44(10):1097–104.

    Article  PubMed  Google Scholar 

  44. Luinenburg DG, de Haan G. MicroRNAs in hematopoietic stem cell aging. Mech Ageing Dev. 2020;189: 111281.

    Article  CAS  PubMed  Google Scholar 

  45. Liu L, Rando TA. Manifestations and mechanisms of stem cell aging. J Cell Biol. 2011;193(2):257–66.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Yang Y, Zhong J, Cui D, Jensen LD. Up-to-date molecular medicine strategies for management of ocular surface neovascularization. Adv Drug Deliv Rev. 2023;201: 115084.

    Article  CAS  PubMed  Google Scholar 

  47. Skeens HM, Brooks BP, Holland EJ. Congenital aniridia variant: minimally abnormal irides with severe limbal stem cell deficiency. Ophthalmology. 2011;118(7):1260–4.

    Article  PubMed  Google Scholar 

  48. Abtahi MA, Beheshtnejad AH, Latifi G, Akbari-Kamrani M, Ghafarian S, Masoomi A, et al. Corneal epithelial thickness mapping: a major review. J Ophthalmol. 2024;2024:6674747.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Abusamak M. Corneal epithelial mapping characteristics in normal eyes using anterior segment spectral domain optical coherence tomography. Transl Vis Sci Technol. 2022;11(3): 6.

    Article  PubMed  PubMed Central  Google Scholar 

  50. Wang X, Dong J, Wu Q. Corneal thickness, epithelial thickness and axial length differences in normal and high myopia. BMC Ophthalmol. 2015;15:49.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Lobo EP, Delic NC, Richardson A, Raviraj V, Halliday GM, Di Girolamo N, et al. Self-organized centripetal movement of corneal epithelium in the absence of external cues. Nat Commun. 2016;7: 12388.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Puri S, Sun M, Mutoji KN, Gesteira TF, Coulson-Thomas VJ. Epithelial cell migration and proliferation patterns during initial wound closure in normal mice and an experimental model of limbal stem cell deficiency. Invest Ophthalmol Vis Sci. 2020;61(10): 27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Porpiglia E, Mai T, Kraft P, Holbrook CA, de Morree A, Gonzalez VD, et al. Elevated CD47 is a hallmark of dysfunctional aged muscle stem cells that can be targeted to augment regeneration. Cell Stem Cell. 2022;29(12):1653-68 e8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Jiang D, et al. Spatiotemporal single-cell analysis elucidates the cellular and molecular dynamics of human cornea aging. HRA007528, 2025. https://ngdc.cncb.ac.cn/gsa-human/browse/HRA007528.

  55. Jiang D, et al. Spatiotemporal single-cell analysis elucidates the cellular and molecular dynamics of human cornea aging. HRA007620, 2025. https://ngdc.cncb.ac.cn/gsa-human/browse/HRA007620.

Download references

Acknowledgements

Not applicable.

Funding

This work was supported by the Key Research and Development Plan of Zhejiang Science and Technology Department (Grant no. 2024 C03206) and the National Key Technologies Research and Development Programme of China (Grant no. 2019YFC0840708).

Author information

Author notes

  1. Dan Jiang, Ke Li and Yining Sun contributed equally to this work.

Authors and Affiliations

  1. National Clinical Research Center for Ocular Diseases, State Key Laboratory of Ophthalmology, Optometry and Visual Science, Eye Hospital, Wenzhou Medical University, Wenzhou, 325027, China

    Dan Jiang, Ke Li, Yining Sun, Zicheng Zhang, Shuang Xie, Xintong Yu, Ruoqi Wang, Ying Feng, Qinxiang Zheng, Yajing Wen, Peter S. Reinach, Meng Zhou & Wei Chen

  2. School of Biomedical Engineering, Wenzhou Medical University, Wenzhou, 325027, China

    Ke Li, Zicheng Zhang & Meng Zhou

  3. Hangzhou Institute of Medicine (HIM), Chinese Academy of Sciences, Hangzhou, 310000, China

    Yuanyuan Du

Authors
  1. Dan Jiang
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Ke Li
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Yining Sun
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Zicheng Zhang
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Shuang Xie
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Xintong Yu
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Ruoqi Wang
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Ying Feng
    View author publications

    You can also search for this author inPubMed Google Scholar

  9. Qinxiang Zheng
    View author publications

    You can also search for this author inPubMed Google Scholar

  10. Yajing Wen
    View author publications

    You can also search for this author inPubMed Google Scholar

  11. Peter S. Reinach
    View author publications

    You can also search for this author inPubMed Google Scholar

  12. Yuanyuan Du
    View author publications

    You can also search for this author inPubMed Google Scholar

  13. Meng Zhou
    View author publications

    You can also search for this author inPubMed Google Scholar

  14. Wei Chen
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

WC, MZ and YYD conceived and designed the study. KL and ZCZ performed the scRNA-seq and spatial transcriptomics data analysis. DJ, YNS, SX, XTY, RQW, YJW, QXZ and YF prepared the samples and performed the experiments. KL, DJ, YNS and MZ drafted the manuscript. MZ, YNS, DJ, YYD and PSR revised the manuscript. MZ, WC and DJ supervised the study. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Yuanyuan Du, Meng Zhou or Wei Chen.

Ethics declarations

Ethics approval and consent to participate

This study was conducted in accordance with the principles of the Declaration of Helsinki. The use of corneal tissues was approved by the Ethical Review Committee of the Affiliated Eye Hospital of Wenzhou Medical University (2022-159-K-124-01). Written informed consent was obtained from the donors'legal representatives prior to tissue collection. The mouse experiments were approved by the Laboratory Animal Ethics Committee of Wenzhou Medical University (YSG24032101).

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

13073_2025_1475_MOESM1_ESM.docx

Additional file 1: Fig. S1. Spatiotemporal transcriptomic and cell type atlas of the human cornea. Fig. S2. Spatiotemporal changes in limbal stem cells during corneal aging. Fig. S3. Spatiotemporal changes in the limbal stem cells of the human cornea. Fig. S4. Dynamic changes in molecular profiles during limbal stem cell aging. Fig. S5. Expression of corneal basal cell markers and cluster-specific markers in corneal epithelial cells. Fig. S6. Dynamic changes in molecular profiles during corneal epithelial cell aging. Fig. S7. Diversity and spatial location of corneal basal cells during aging. Fig. S8. Changes in the cellular neighbourhoodsof the corneal epithelium with aging. Fig. S9. Cellular interactions and THBS1-associated gene expression in corneal basal cells. Fig. S10. Dynamic changes in molecular profiles during corneal endothelial aging

13073_2025_1475_MOESM2_ESM.xlsx

Additional file 2: Table. S1: Lists of differentially expressed genes in various cell types. Table. S2: Differentially expressed genes of different anatomical regions in spatial transcriptomics data. Table. S3: Differentially expressed genes of different LSC subsets. Table. S4: Differentially expressed genes of different corneal basal subsets. Table. S5: Differentially expressed genes of different stromal cell subsets

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Jiang, D., Li, K., Sun, Y. et al. Spatiotemporal single-cell analysis elucidates the cellular and molecular dynamics of human cornea aging. Genome Med 17, 56 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13073-025-01475-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13073-025-01475-z

Keywords

Genome Medicine

ISSN: 1756-994X

Contact us