In situ spatial transcriptomics reveals novel markers of the limbal stem cell niche and ocular surface epithelia

TL;DR

This study uses spatial transcriptomics to identify new markers for limbal stem cells and neutrophils in the eye, which could improve treatments for corneal stem cell deficiency.

Contribution

The study identifies Krt16 and Nkiras1 as novel markers for limbal stem cells and neutrophils using spatial transcriptomics.

Findings

Krt16 is dynamically expressed in the developing limbus and during corneal injury, marking functional stem cells.

Nkiras1 identifies a population of limbal neutrophils.

Limbal span forms in mice postnatally after eyelid opening.

Abstract

The mammalian cornea is endowed with stem cells (SCs) that have lifelong regenerative activity. The niche for these cells is the limbus, and damage to it or its SCs results in limbal stem cell deficiency (LSCD). Despite the numerous studies that employ single-cell RNA sequencing, the identity of these cells remains an enigma principally because their spatial positioning is lost upon dissociation. These adversities were avoided via on-tissue spatial transcriptomics where Krt16 and Nkiras1 were differentially expressed. Krt16 was dynamically expressed in the developing limbus, correlated with slow-cycling label-retaining limbal epithelial SCs and was induced during corneal injury, observations consistent with marking functional SCs. Additionally, we established Nkiras1 as a novel maker of limbal neutrophils. Because current gold-standard treatments for LSCD include SC transplantation, our data will inform future studies in delivering a more reliable standard therapy that incorporates an identifiable SC population to improve clinical outcomes.

•Spatial transcriptomics on ocular surface tissue identified K16 as a limbal stem cell marker•K16+ epithelial cells partake in corneolimbal development and wound healing•NKIRAS1+ cells identify a population of limbal neutrophils•Limbal span forms in mice postnatally after eyelid open

Spatial transcriptomics on ocular surface tissue identified K16 as a limbal stem cell marker

K16+ epithelial cells partake in corneolimbal development and wound healing

NKIRAS1+ cells identify a population of limbal neutrophils

Limbal span forms in mice postnatally after eyelid open

In this article, Di Girolamo and colleagues report on the tissue-specific spatial transcriptomic profiles of healthy corneal, stem cell-containing limbal and conjunctival epithelia using the mouse as a model. Unique genes and networks are identified and validated that may improve cell isolation for therapeutic use and assist the ongoing search for reliable stem cell biomarkers to enhance clinical outcomes.

Figures

Differential expression analysis of genes across selected ocular surface regions(A) (i) Representative image of a hematoxylin-stained cross-section through the mouse ocular surface. Scale bar represents 100 μm. (ii) Representative image of K8/K12 double-stained tissue cross-section taken at 20x magnification using a scanning confocal microscope. Scale bars represent 50 μm. (iii) Representative image of ocular surface tissue cross-section immunostained with K5 and DAPI for *in situ* spatial profiling captured through the GeoMx DSP software. Hatched rectangles represent ROIs chosen for transcrip

Signaling patterns and communication networks across the ocular surface(A) Scatterplot displays the relationship between incoming and outgoing interactional strength across different regions of the ocular surface. Dot size represents relative contribution.(B) Dot plot depicts outgoing and incoming signaling patterns between tissue regions. Color intensity represents signaling strength; darker colors indicate stronger signals; lighter colors indicate weaker signals.(C) Circular plots indicate specific gene regulatory networks enriched across the ocular regions. Connecting lines between circles

Validating selected DEGs via molecular and immunohistochemical methods(A) Representative graphs summarize gene expression for corneal *Krt12* and *Slurp1* (i & ii); limbal *Krt16*, *Krt17*, *Apoe*, and *Nkiras1* (iii–vi); and conjunctival *Krt8* and *Alox12e* (vii and viii) genes by qPCR. Violin plots (first column) represent average expression level of targets identified via the GeoMx spatial transcriptomics assay (*n* = 3 in triplicate). Second and third columns represent normalized mRNA fold-change between two regions. qPCR data were obtained from 8 replicates (pooled from *n* = 16 mice and

K16 identifies a subset of limbal epithelia with important physiological roles(A) Representative images of K16-stained epithelia on whole flat-mounted ocular surface tissue, captured at 63x magnification (*n* = 4). Maximum intensity projection images of corneal, limbal, and conjunctival epithelia are displayed in (a–c). Scale bars represent 20 μm. Analysis of K16 fluorescence intensity (*n* = 4) is shown in (d). Enlarged images of K16 expression in segmented basal/suprabasal and superficial layers are shown in (e and f). Scale bars represent 20 μm.(B) Representative images of K16 staining on t

NKIRAS1 as a novel marker of neutrophils and their functional significance in physiological processes(A) Double immunostaining for K8 and NKIRAS1 on whole flat-mounted tissue (*n* = 4) captured at 20x magnification. Hatched white line represents the limbal boundary to the cornea. Scale bars represent 50 μm.(B) Representative images of NKIRAS1^+^ cells on whole flat-mounted ocular surface tissue captured at 63x magnification, showing multilobular nuclei associated with NKIRAS1^+^ cells (*n* = 5). Scale bars represent 20 μm.(C) Representative images of NKIRAS1/LY6G double immunostaining on whole

Full text

Introduction

Optimal visual perception is achieved by the maintenance of a healthy ocular surface, which includes the cornea, limbus, and conjunctiva. The narrow limbal zone, which intersects the cornea and conjunctiva, is adorned with limbal epithelial stem cells (LESCs) (Cotsarelis et al., 1989; Pellegrini et al., 1999). These cells have self-renewal and tissue regenerative capacity by replenishing the stem cell (SC) pool and replacing aged or damaged epithelia exfoliated from superficial layers (Davanger and Evensen, 1971; Di Girolamo, 2015), while acting as a barrier to conjunctival incursion (Di Girolamo and Park, 2023; Dua and Azuara-Blanco, 2000). Depletion of LESCs or damage to the limbus triggers a vision-threatening condition known as limbal stem cell deficiency (LSCD), whereby corneal regeneration is disrupted and pathological conjunctival cells invade into corneal terrain (Deng et al., 2019).

Currently, SC transplantation in various formats is the treatment of choice to restore sight in LSCD. However, despite significant therapeutic advances (Ghareeb et al., 2020), patients continue to suffer poor long-term outcomes, including low visual acuity and persistent neovascularization (Behaegel et al., 2019; Rama et al., 2010; Zhao and Ma, 2015), ultimately reducing quality of life. The limbal frontier with the cornea and conjunctiva is well characterized by structural landmarks (Nureen et al., 2024; Nureen and Di Girolamo, 2024); however, epithelia within this zone are poorly defined. This knowledge gap has stalled the development of efficacious strategies for ocular surface reconstruction because biomarkers that identify LESCs remain concealed. This insight is key to facilitating their isolation and detailing how they restore tissue architecture and sight upon transplantation. The expression of structural proteins such as keratins (K) including K3/K12 for corneal (Liu et al., 1994; Rodrigues et al., 1987) and K4/K13 or K8 for conjunctival (Krenzer and Freddo, 1997; Ramirez-Miranda et al., 2011) epithelia is well established. Although K14, K15, and K19 are regarded markers of limbal epithelia, they also display coincidental expression in conjunctival cells (Chen et al., 2004; Kao, 2020; Krenzer and Freddo, 1997; Yoshida et al., 2006). Putative markers of LESCs, including ΔNP63, ABCB5, LRIG1, K14, K15, GPHA2, IFITM3, GAS1, and HES1 identified by in situ hybridization (Di Iorio et al., 2005), quantitative lineage tracing (Altshuler et al., 2021), in situ photo-labeling (Farrelly et al., 2021), and transgenic animal studies (Ksander et al., 2014; Nakamura et al., 2008; Richardson et al., 2025), can also be found in neighboring epithelia (Nureen et al., 2024), further complicating their identification. More recently, transcriptomic profiling of the human and mouse ocular surface via high-throughput single-cell RNA sequencing (scRNA-seq) elucidated differentially expressed genes (DEGs) in the limbus compared to adjacent locations (Collin et al., 2021; Dou et al., 2021; Kaplan et al., 2019; Li et al., 2021; Ligocki et al., 2021; Lin et al., 2023). This strategy is popular for biomarker discovery and characterizing cell subpopulations. However, there are caveats to consider, including sample heterogeneity, mode of cornea procurement and storage, exposure to chemical and physical stressors, loss of spatial resolution, and sequencing depth. These factors can alter the original transcriptional landscape of the tissue and impose variability between studies, rendering datasets incompatible and unreliable (Arts et al., 2023).

Herein, we exploited an on-tissue spatial transcriptomics approach to overcome limitations of prior methodologies. This platform enabled accurate identification of epithelial cells within the corneal, limbal, and conjunctival precincts based on morphology and keratin expression, while preserving the tissue for gene expression profiling at the whole transcriptome level. Our precision mapping of the limbal expanse with faithful markers (Nureen et al., 2024; Nureen and Di Girolamo, 2024) lends credence to performing the first spatial sequencing of ocular surface cells without perturbing their position. These findings will broaden our understanding of LESC regulation, function, and interplay with neighboring cells and their niche under steady-state, wound-healing, and disease-baring conditions.

Results

Differential keratin protein expression defines ocular surface regions

To accurately delineate each epithelial division for spatial transcriptomic profiling, serial ocular surface tissue sections from mouse eyes were stained with hematoxylin (Figure 1A(i)) and double immunostained with K8 and K12 (Figure 1A(ii)). The small zone harboring a dual-layered epithelium, devoid of K8/K12 immunoreactivity, was confidently deemed the murine limbus (Nureen et al., 2024; Nureen and Di Girolamo, 2024). This survey facilitated the selection of four distinct regions of interest (ROIs), including a threshold of 20–25 nuclei counts for corneal (K12^+^), limbal (K8^−^/K12^−^), and conjunctival (K8^+^) epithelia (Figure 1A(ii)). A “transitional” zone with the same threshold nuclei count was incorporated to capture epithelial cells arising from the limbal neighborhood, likely to include progenitor/early transient amplifying cells (TACs) (Figure 1A(ii) and (iii)). K5 was utilized as a fluorescent morphology marker to detect all epithelial cells within individual tissue sections (Figure 1A(iii)).Figure 1. Differential expression analysis of genes across selected ocular surface regions(A) (i) Representative image of a hematoxylin-stained cross-section through the mouse ocular surface. Scale bar represents 100 μm. (ii) Representative image of K8/K12 double-stained tissue cross-section taken at 20x magnification using a scanning confocal microscope. Scale bars represent 50 μm. (iii) Representative image of ocular surface tissue cross-section immunostained with K5 and DAPI for in situ spatial profiling captured through the GeoMx DSP software. Hatched rectangles represent ROIs chosen for transcriptomic profiling. Scale bar represents 100 μm.(B) Heatmap shows expression pattern of marker genes for each ROI across the four ocular regions (color-coded vertical columns). Red indicates high and blue indicates low expression.(C) Volcano plots show markers from DEA of genes for limbus (i), cornea (ii), conjunctiva (iii), and transitional region (iv) compared to all other regions. Hatched boxes are magnified views of key DEGs. x axis shows log2 fold-change, and y axis shows −log10 p value. DESeq2 was used to determine DEGs, and significance was assessed after correcting multiple hypotheses testing with the Benjamini-Hochberg procedure. All data and images were obtained from n = 3 independent eyes of different mice in triplicate tissue sections.

Region-specific transcriptional signatures

Next, we employed the GeoMx Digital Spatial Profiling (DSP) platform using a whole-transcriptome gene panel to ascertain gene expression in the selected tissue divisions. After quality control and data processing, 34 ROIs were retained for downstream analysis. The limbus was the focus of our investigation because it is purported to harbor SCs. Here, we identified eight upregulated genes in the limbus vs. all other regions within the ocular surface (ordered by magnitude of average log2 fold-change), Krt16, Krt17, Csnk2a2, Trim46, Mccc2, Crtam, Nkiras1, and Apoe, and one downregulated gene Gapdh (Figures 1B and 1C(i)). Spatial maps of selected DEGs are displayed in Figure S1. Krt16 is noteworthy due to its role as a stress- or damage-associated keratin that regulates innate immunity in response to epidermal barrier breach (Lessard et al., 2013), and it is also a marker of prostate epithelial SCs (Hu et al., 2021). The other keratin identified as a limbal-specific DEG was Krt17 (Figures 1B and 1C(i)). These data are congruent with its expression in undifferentiated cells (Bonnet et al., 2021) and recognition as a marker of putative LESCs (Lin et al., 2023; Wang et al., 2025). A novel limbal DEG that caught our attention was Nkiras1, a negative regulator of nuclear factor (NF)-κB signaling, preventing NF-κB degradation, thereby modulating inflammatory networks (Gerashchenko et al., 2010). In the current setting, it may serve multiple roles including protection against inflammation-induced tissue damage and suppression of tumorigenesis (Postler et al., 2023) that can arise in SC populations. Apoe (apolipoprotein E) plays a role in lipoprotein metabolism and neuroprotection, and its loss of function increases susceptibility to oxidative stress and neurodegeneration (Shea et al., 2002). In harmony with our data, Apoe is considered a putative marker of LESCs (Li et al., 2024; Lin et al., 2023). Gapdh, a classic glycolytic enzyme, was the only downregulated DEG in limbus compared to other regions. This likely reflects the metabolic quiescence of SCs within the limbal compartment (Altshuler et al., 2021; Kulkarni et al., 2010).

In the corneal epithelium, 73 upregulated and 12 downregulated genes were identified (Figures 1B and 1C(ii)) when compared to the other three regions. The most significantly upregulated gene was Slurp1 (secreted Ly-6/uPAR-related protein-1; avg_log2FC +2.5) (Figures 1B and 1C(ii)). Slurp1 is a secreted immunomodulatory protein involved in epithelial homeostasis, and its overexpression in human corneal epithelia increases p15/CDKN2B, a cyclin-dependent kinase inhibitor, suggesting that it is an anti-proliferative factor (Swamynathan et al., 2022) that regulates angiogenic inflammation (Swamynathan et al., 2024). Other upregulated DEGs in the corneal epithelium included Krt12, Lypd2, Erich5, Tkt, and Aqp5. Consistent with our results, many of these genes have been detected in the cornea (Ligocki et al., 2021; Lin et al., 2023; Sax et al., 1996), thus validating our in situ transcriptomic assay. Krt12 is highly specific to differentiated corneal epithelia (Nureen et al., 2024; Nureen and Di Girolamo, 2024; Wolosin, 2024). As such, it effectively discriminated corneal from conjunctival epithelia, with an avg_log2FC of −5.53.

In the conjunctival epithelium, 35 upregulated and 25 downregulated genes were identified (Figures 1B and 1C(iii)). The most highly expressed DEGs were Alox12e, Alox15, Krt4, Krt7, Krt8, Krt19, Cbr2, and Cyp2f2 (Figure 1C(iii)). Both Alox12e and Alox15 belong to the lipoxygenase family of enzymes involved in arachidonic acid metabolism and production of bioactive lipid mediators. Alox12e catalyzes 12-hydroxyeicosatetraenoic acid, a molecule that regulates inflammation, leukocyte recruitment, and vasoconstriction (Kulkarni et al., 2021). We also detected Krt8 as a highly expressed conjunctival gene, validating the robustness and reliability of our spatial GeoMx transcriptomics data. Notably, cornea-specific Slurp1, Lypd2, Tkt, and Aqp5 were significantly downregulated in this domain (Figure 1C(iii)).

The transitional zone consisting of the peripheral cornea and limbus exhibited fewer unique DEGs (Figure 1C(iv)). This is because the ROIs selected likely captured segments of cornea and limbus, as we attempted to characterize a region that contains TACs. Olfr716, Rfng, Mocs3, and Pard6b were among the induced DEGs in this region. Pard6b (partitioning defective 6 homologue beta) is a regulator of cell polarity (Alarcon, 2010) and cell-cycle restrictor (Marques et al., 2016), functions that resonate with epithelia proximal to the limbus. The three downregulated DEGs (Alox12e, Krt7, Krt19) are conjunctival specific (Figure 1C(iii)) and expected to be suppressed in this district (Lin et al., 2023).

Cell-to-cell communication networks within the ocular surface

Communication networks across the four ocular surface landscapes focused on cell-to-cell contact, extracellular matrix-receptor interactions (Figure 2A), strength of the interactions, and secreted signaling pathways. This analysis identified 39 distinct gene regulatory networks governing ocular surface homeostasis and region-specific functions (Figure 2B). The conjunctiva emerged as the core for both incoming and outgoing signaling interactions (Figure 2A). Specifically, it was deemed the primary hub for immune interactions and antigen processing, with significant involvement in neutrophil-mediated inflammatory responses through CEACAM1 signaling (Figure 2C). This conjunctiva-centric signaling resonates with its critical role in mucosal immunity and immune surveillance (de Paiva et al., 2022). The limbus plays a central role in orchestrating immune responses, particularly through CCL chemokines and interleukin (IL)-2 signaling (Figure 2C). It engages with the transitional region via CCL chemokine signaling, which is essential for neutrophil chemotaxis into inflammatory sites across the ocular surface (Metzemaekers et al., 2020; Wallace et al., 2004). In parallel, the limbus communicates with the conjunctiva through IL-2 signaling, a key component of “stem-microenvironmental networking.” IL-2 bridges innate and adaptive response by modulating the activity of immune cells, including T cells and natural killer cells.Figure 2. Signaling patterns and communication networks across the ocular surface(A) Scatterplot displays the relationship between incoming and outgoing interactional strength across different regions of the ocular surface. Dot size represents relative contribution.(B) Dot plot depicts outgoing and incoming signaling patterns between tissue regions. Color intensity represents signaling strength; darker colors indicate stronger signals; lighter colors indicate weaker signals.(C) Circular plots indicate specific gene regulatory networks enriched across the ocular regions. Connecting lines between circles indicate signal interactions between regions, and colors indicate the provenance of the interaction. Line thickness represents interactional strength.All data were obtained from n = 3 independent eyes of different mice in triplicate tissue sections.

The transitional zone also exhibits involvement in immune surveillance through major histocompatibility complex (MHC) class I expression (Figure 2C). MHC class I presents antigens to cytotoxic T cells, playing a crucial role in identifying and eliminating diseased cells. MHC class I expression has been reported in various ocular surface cells (Hamrah et al., 2002), indicating this region’s role in immune vigilance. Collectively, these findings underscore the complex molecular interactions that maintain ocular surface integrity and highlight region-specific functions essential for effective immune protection and response.

Pathway enrichment within the ocular surface

Pathway enrichment analysis revealed a significant role for the corneal and transitional zones in xenobiotic metabolism (Figure S2), a pivotal activity for clearing environmental and microbial toxins and carcinogens. Within this pathway, the expression of several genes was increased in the corneal compared to transitional compartment. Notably, prostaglandin reductase-1 (Ptgr1) and inhibitor of differentiation-2 (Id2) emerged as corneal biomarkers through differential expression analysis (DEA) (Figure 1C(ii)). Ptgr1 is modulated in damaged corneas (Wang et al., 2024b), suggesting a role in maintaining tissue integrity. Id2 is differentially expressed by corneal epithelial cells and corneal fibroblasts (Mohan et al., 2016). Additionally, the presence of limbal-associated Mccc2 and Apoe suggests that a collaborative interaction between the cornea and limbus exists to facilitate detoxification (Figure S2A).

The transitional zone is enriched with genes downregulated by KRAS activation, including Krt4 and Krt15 (Figure S2B), which are established markers of conjunctival and limbal epithelia (Altshuler et al., 2021; Harun et al., 2013; Nasser et al., 2018). KRAS plays a central role in signal transduction, initiating RAF/MEK/ERK and PI3K/AKT signaling cascades. Dysregulation in this pathway is associated with tumorigenesis and metastases (Huang et al., 2021; Lam et al., 2022), which are unwanted processes in a SC-harboring location. The unique gene expression profile of the transitional region, influenced by its proximity to both limbus and cornea, suggests a specialized role in modulating cellular proliferation and maintaining tissue homeostasis.

Validating DEGs by qPCR

Selected DEGs detected in each tissue division of the murine ocular surface were validated by reverse-transcription-quantitative polymerase chain reaction (RT-qPCR) (Figure 3A). Because accurate tissue procurement from the transitional region was technically challenging, it was excluded from further investigations. Increased expression for Krt16 (7.74-fold), Krt17 (72.73-fold), Apoe (7.27-fold), Trim46 (4.23-fold), and Crtam (16.67-fold) was observed in the limbus compared to cornea (Figures 3A(iii–v) and S3). Expression of other limbal DEGs, such as Mccc2, Nkiras1, and Csnk2a2, also increased by 21.38-, 4.67-, and 15.54-fold in the limbus compared to cornea, respectively (Figure S3). Nkiras1 mRNA was significantly increased (27.99-fold) (Figure 3(vi)) when comparing limbus to conjunctiva. The remaining DEGs, Krt16, Krt17, Mccc2, Crtam, and Csnk2a2, were increased by 12.39-, 66.1-, 9.24-, 2.90-, and 12.39-fold in the limbus compared to conjunctiva (Figures 3A(iii and iv) and S3). Unexpectedly, Apoe and Trim46 mRNAs were increased by 1.66- and 2.07-fold in the conjunctiva vs. limbus, respectively (Figures 3A(v) and S3). Although Gapdh was identified as a downregulated limbal DEG, this was not the case by qPCR (Figure S3). This discrepancy may have arisen due to the difficulty in dissecting tissues without contamination from adjacent sites.Figure 3. Validating selected DEGs via molecular and immunohistochemical methods(A) Representative graphs summarize gene expression for corneal Krt12 and Slurp1 (i & ii); limbal Krt16, Krt17, Apoe, and Nkiras1 (iii–vi); and conjunctival Krt8 and Alox12e (vii and viii) genes by qPCR. Violin plots (first column) represent average expression level of targets identified via the GeoMx spatial transcriptomics assay (n = 3 in triplicate). Second and third columns represent normalized mRNA fold-change between two regions. qPCR data were obtained from 8 replicates (pooled from n = 16 mice and 32 eyes) and analyzed using unpaired Student’s t test. Error bars represent ±SEM. ^∗^p < 0.05, ^∗∗^p < 0.01; ^∗∗∗∗^p < 0.0001.(B) Representative images of RNA expression for selected DEGs in cornea, limbus, and conjunctiva detected by RNAscope (n = 3). Images were captured at 63x magnification using a scanning confocal microscope. Negative control includes samples without target probes, which developed no hybridization signal. DAPI (blue) counterstains the cell nuclei. Hatched white line partitions the epithelium from stroma. Scale bars represent 20 μm.(C) Representative images of selected marker staining by immunofluorescence on whole flat-mounted ocular surface spanning cornea to conjunctiva (n = 5 per marker). Hoechst (blue) counterstained the cell nuclei. Images were captured with 20x objective lens. Hatched white lines demarcate limbal boundaries with cornea (inner) and conjunctiva (outer), based on basal cell nuclei morphology and staining for K8 and K12. Scale bars represent 50 μm.

Additionally, two corneal and two conjunctival genes were chosen for qPCR-based validation. As expected, corneal Krt12 and Slurp1 transcripts were significantly increased by 920.73- and 494.31-fold in the cornea compared to conjunctiva (Figure 3A(i and ii)), corroborating previous scRNA-seq data (Lin et al., 2023). Similarly, Krt8 and Alox12e were amplified in the conjunctival compartment by 75,636.30- and 50,076.10-fold compared to cornea, and by 23.78- and 12.78-fold compared to limbus, respectively (Figure 3A(vii and viii)).

Localization of selected DEGs by RNAscope and immunofluorescence

RNAscope revealed specific punctate staining for Krt16 and Apoe reflecting RNA expression in limbal epithelia (Figure 3B), corroborating to the findings from the GeoMx spatial assay (Figure 1C). Alox12e was also validated as a novel and specific marker of conjunctival epithelia (Figure 3B). Krt8 and Krt12 RNA transcripts were confined to corneal and conjunctival epithelia, respectively (Figure 3B), which is consistent with the pattern of K8 and K12 protein expressions, observed by immunofluorescence (Figure 3C). K16 protein was detected in discrete clusters of limbal epithelia (Figure 3C). Nkiras1 expression was investigated due to its novelty on the ocular surface and established role as a tumor suppressor (Postler et al., 2023). It was minimally expressed by limbal epithelia (Figure S4) but predominated in cells within the murine limbal and conjunctival stroma (Figure 3C).

Next, K16 protein was elaborated on whole flat-mounted ocular surface tissue where it was detected in the limbus, absent from the stroma, negligible in cornea, and low in conjunctiva (Figure 4A(a–c)). Fluorescence intensity measurements confirmed higher K16 expression in the limbal compared to corneal (p = 0.046) and conjunctival (p = 0.009) epithelia (Figure 4A(d)). Notably, K16 staining spanned basal-to-superficial limbal epithelia (Figures 4A(e and f) and S5), an observation supported by recent scRNA-seq data (Lin et al., 2023; Wang et al., 2025), suggesting that these cells possess both proliferative and regenerative potential. K16 was highly expressed in human basal limbal epithelia but minimally stained corneal epithelial cells (Figure 4B), thereby corroborating data accrued from mice. However, a small number of stromal cells expressed K16 on human tissues (Figure 4B). The prevailing notion is that keratins are the hallmark of epithelial cells. Notably, these proteins can also be expressed by T cells (He et al., 2025) and mesenchymal progenitor cells (Ferretti et al., 1989) both of which exist in the limbal stroma. Curiously, macrophages (Lu et al., 2024), melanocytes (Bhawan et al., 2005), and cells undergoing epithelial-to-mesenchymal transition (Elazezy et al., 2021) can also express K16. Given our human corneas were resected from donor globes and stored in hypothermic solution with a postmortem delay of <64 h, it is not surprising that a stress-related keratin was detected in the stromal compartment (Figure 4B). However, the identity of these cells remains unknown.Figure 4K16 identifies a subset of limbal epithelia with important physiological roles(A) Representative images of K16-stained epithelia on whole flat-mounted ocular surface tissue, captured at 63x magnification (n = 4). Maximum intensity projection images of corneal, limbal, and conjunctival epithelia are displayed in (a–c). Scale bars represent 20 μm. Analysis of K16 fluorescence intensity (n = 4) is shown in (d). Enlarged images of K16 expression in segmented basal/suprabasal and superficial layers are shown in (e and f). Scale bars represent 20 μm.(B) Representative images of K16 staining on the human cornea and limbus (n = 4). Scale bars represent 50 μm.(C) Representative images of BrdU (green) immunofluorescence staining on whole-mount tissue captured at 40x magnification (n = 4). Hatched white lines represent limbal boundaries to the cornea and conjunctiva. Scale bar represents 20 μm.(D) Line plot of K16 and BrdU^+^ cell distribution across quadrants (n = 4/group).(E) Comparison of the number of K16^+^ cells and BrdU^+^ nuclei (n = 4/group).(F) Correlation plot of average count of K16^+^ and BrdU^+^ cell count (n = 4/group) using linear regression model.(G–L) (G) K16 expression in whole cornea and different regions of wounded vs. unwounded cornea (n = 4/group) at 20x magnification. Scale bars represent 500 μm in first column and 30 μm in other columns. K16^+^ areas in each region are shown in (H), and total percentage areas are shown in (I and J). Representative images of K16 expression in different regions within whole corneas from 1-week-old mice (n = 4) vs. 8-week-old mice (n = 5) at 20x magnification. Scale bars represent 500 μm in first column and 50 μm in other columns. Mean K16^+^ area in each region is shown in (K), and total percentage area is shown in (L); Hoechst (blue) counterstains cell nuclei.All data were analyzed using Student’s t test. Bars represent mean ± SEM. Co, cornea; Conj, conjunctiva; ns, not significant; ^∗^p < 0.05, ^∗∗^p < 0.01, ^∗∗∗^p < 0.001.

Correlation between K16+ epithelia and BrdU+ label-retaining cells

Because label-retaining cells (LRCs) likely represent SCs (Figure 4C) (Nureen et al., 2024; Zhao et al., 2009), we explored the possibility that K16 marks these cells. Because our anti-K16 antibody does not react in acidic solutions (necessary for bromodeoxyuridine [BrdU] detection), the inability to visualize the co-expression of BrdU/K16 double-positive cells is a limitation of our study. Future label-retaining studies will consider using a less detrimental reagent such as 5-ethynyl-2'-deoxyuridine (EdU) to detect slow-cycling SCs. Nonetheless, we determined whether a correlation exists between cells expressing either marker. As expected, no significant difference between the average number of K16^+^ and BrdU^+^ cells (p = 0.359) (Figure 4D) and ocular surface quadrant location (p > 0.05) was detected (Figure 4E) with a correlation coefficient of R^2^ = 0.869 between these two cell groups (Figure 4F). Interestingly, the K16^+^ cell counts appeared slightly higher in the inferior temporal (IT) quadrant (Figure 4D), although the functional relevance of this observation remains unresolved.

K16 expression during corneal wound healing and post-natal development

Next, we interrogated whether K16^+^ limbal epithelial cells contribute to corneal wound healing. In this scenario, K16 was significantly upregulated in the injured central (p < 0.001) and peripheral cornea (p = 0.019) and limbus (p = 0.014) compared to uninjured contralateral eyes, 24 h post-induction of a 2-mm central corneal epithelial debridement (Figures 4G and 4H). Accordingly, the total area of K16^+^ staining in wounded eyes was significantly increased compared to controls (p = 0.002) (Figure 4I), a result confirmed by qPCR (Figure S7). Re-epithelialization at this early time point featured elongated cells that were arranged as linear tracks, consistent with clones spawning from the limbus then moving and converging to form wound-closure lines (Figure 4G) (Dua and Forrester, 1987).

To confirm the hyperproliferative capacity of K16^+^ cells (Schermer et al., 1989), we next examined their distribution in neonatal (1 week old) mouse corneas, a paradigm that reflects a state of heightened growth compared to 8-week-old mice when this activity is subdued (Richardson et al., 2017). Notably, in new-born mice, the boundary separating the K8^+^ conjunctival from K12^+^ corneal zone is blurred compared to the obvious partition that exists in adulthood (Figure S8). Rather, epithelia on either side of this divide intermingle and jostle for position (white arrow; Figure S8). As such, K16 expression is heightened in the central (p = 0.416) and peripheral cornea (p = 0.012) and limbal region (p = 0.003) in 1-week- compared to 8-week-old mice (Figures 4J and 4K). The total area of K16^+^ staining was significantly higher at 1 week compared to 8 weeks (Figure 4L), suggesting that prior to eyelid opening, K16^+^ progenitor cells segregate into the niche to secure a position and impose their epithelial barrier function.

Nkiras1 marks neutrophils in the limbal stem cell niche

Nkiras1 protein was detected in the limbal-conjunctival boundary within the stroma, denoted by the conjunctival marker, K8 (Figure 5A). The multilobulated morphology of NKIRAS1^+^ cells (Figure 5B) and their stromal location are features indicative of neutrophils. This assumption was confirmed by the co-expression of LY6G (Calcagno et al., 2021) and NKIRAS1 (Figure 5C). Notably, there were more NKIRAS1^+^/LY6G^+^ compared to NKIRAS1 (p < 0.001) or LY6G (p < 0.001) single-positive cells (Figure 5D). Accordingly, NKIRAS1^+^/LY6G^+^cells (88.50% ± 2.04%) were numerous compared to single-positive LY6G^+^ (11.50% ± 2.04%) or NKIRAS1^+^ (0%) cells (Figure 5E). Moreover, NKIRAS1^+^/LY6G^+^ neutrophils were irregular shaped and appeared larger within the limbus compared to conjunctiva (Figure S9). The differential immunostaining and morphology suggest that multiple neutrophil subpopulations exist in this tissue as witnessed by other researchers (Ng et al., 2025; Qi et al., 2021). Notably, NKIRAS1^+^/LY6G^+^ neutrophils were situated near blood vessels at the limbal-conjunctival junction (Figure 5F).Figure 5NKIRAS1 as a novel marker of neutrophils and their functional significance in physiological processes(A) Double immunostaining for K8 and NKIRAS1 on whole flat-mounted tissue (n = 4) captured at 20x magnification. Hatched white line represents the limbal boundary to the cornea. Scale bars represent 50 μm.(B) Representative images of NKIRAS1^+^ cells on whole flat-mounted ocular surface tissue captured at 63x magnification, showing multilobular nuclei associated with NKIRAS1^+^ cells (n = 5). Scale bars represent 20 μm.(C) Representative images of NKIRAS1/LY6G double immunostaining on whole flat-mounted tissue (n = 6) at 40x magnification. White arrows point to LY6G-single-positive cells. Scale bars represent 20 μm.(D) Comparing mean cell density of single- and double-positive NKIRAS1/LY6G-expressing cells (n = 6).(E) Comparing mean percentage of single- and double-positive NKIRAS1/LY6G-expressing cells (n = 6). Data were analyzed by one-way ANOVA.(F) Triple immunostaining of NKIRAS1, LY6G, and phalloidin (for blood vessels) on whole-mount corneas (n = 3). For clarity, phalloidin staining is depicted by the white trace. Scale bars, 50 μm.(G) Nkiras1 expression in whole cornea and different regions of wounded vs. unwounded cornea (n = 4/group) at 20x magnification. Scale bars represent 500 μm in first column and 50 μm in other columns. Average number of NKIRAS1^+^ cells in different regions is shown in (H), and total number is shown in (I). Hoechst (blue) counterstains cell nuclei.Data were analyzed using paired Student’s t test. Bars represent mean ± SEM. ns, not significant; ^∗^p < 0.05, ^∗∗∗^p < 0.001, ^∗∗∗∗^p < 0.0001.

Nkiras1 expression after corneal wounding

Neutrophil recruitment is a key component of acute inflammation that facilitates corneal wound healing (Li et al., 2006; Sumioka et al., 2021). As such, NKIRAS1^+^ neutrophils were elevated in the limbus of wounded compared to uninjured controls (p < 0.001) but not in the central (p > 0.99) or peripheral cornea (p = 0.415) (Figures 5G and 5H). This is not surprising, given that monitoring was conducted over 24 h, a time frame when the wound is sealed and neutrophils begin to subside (Chen et al., 2024). Overall, the number of NKIRAS1^+^ neutrophils was higher in wounded compared to controls (p = 0.020) (Figure 5I), confirming they play a functional role in innate immunity and inflammation.

Discussion

To our knowledge, on-tissue spatial transcriptomic profiling of the mammalian ocular surface with a granular focus on the epithelium has not been conducted. One recent study employed a similar spatial approach with probe-based hybridization technology on human corneas; however, tissue boundaries were not specified and postmortem delay may have influenced the data, especially given the markers uncovered for LESCs were non-specific (Jiang et al., 2025). Most scRNA-seq studies conducted to date employ strategies that expose cells and tissues to enzymes that modify gene expression, phenotype, and function and are performed from poorly defined precincts. Therefore, our goal was to accurately define the transcriptional signature within the confines of the cornea, limbus, and conjunctiva, targeting the limbal compartment, primarily because it is deemed the repository for LESCs, and hence the cornea’s regenerative hub (Cotsarelis et al., 1989; Davanger and Evensen, 1971; Di Girolamo, 2015). Eight upregulated and one downregulated gene was uncovered in this division compared to cornea, transition zone, and conjunctiva (Figure 1). Of these, two genes (Krt16 and Nkiras1) emerged as noteworthy and were prioritized for further investigation. Several new corneal and conjunctival markers were also identified (Figure 1C), thus expanding the transcriptomic atlas for these tissues.

A key finding was the localization of K16 in the SC-containing limbal compartment, specifically within distinct clusters reminiscent of clones (Figures 4A and S5) that were associated with BrdU^+^ LRCs (Figures 4D–F) and its presence in human samples. The spatiotemporal evolution of K16^+^ cells during post-natal life (i.e., a time when the limbal niche undergoes dynamic structural modifications) (Figure S8), and the increased presence of these cells during wound healing (Figure 4G) add credence to its role as a stress-induced tissue-regenerative keratin (Lessard et al., 2013). Notably, K16 is involved in wound repair and re-epithelialization in inflammatory skin disease (Mazzalupo et al., 2003; Patel et al., 2006; Zhang et al., 2019), and Krt16-knockout mice develop skin and oral lesions soon after birth (Lessard and Coulombe, 2012). In this setting, and upon exposing the epidermis to chemical or mechanical irritation, keratinocytes lacking K16 are vulnerable to immune-mediated reactions because critical inflammatory checkpoints are erased or dysfunctional (Lessard et al., 2013). Krt16-null mice have no reported ocular phenotype (Lessard and Coulombe, 2012; Zieman and Coulombe, 2018), possibly due to redundancy imposed by related keratin species such as K17, which we (Figure 3A(iv)) and others (Lin et al., 2023; Wang et al., 2025) discovered in the limbal compartment. Notably, Krt17 deletion results in enrichment of biological processes responsible for pro-inflammatory cytokine response and the stimulation of neutrophil chemotaxis (Wang et al., 2025). Contextually, this resonates with Krt17 being a corneal-immunosuppressive gene. Likewise, our data on Krt16 suggest its expression in limbal progenitor cells may be a key requirement for maintaining an immune-silenced microenvironment in a remit that is constantly exposed to antigenic stimuli. K16 and K17 have also been ascribed unexpected non-mechanical functions. Through their interactions with intracellular organelles (such as mitochondria), they regulate redox homeostasis, thus protecting against oxidative stress (Steen et al., 2020). The strategic positioning of K16^+^ cells within the limbus under steady state (Figures 3C and 4A) suggests a population of SCs is being marked. This is a reasonable assumption since K16 is expressed by prostate stem/progenitor cells (Hu et al., 2021) and marks basal skin keratinocytes (Paladini et al., 1996) and its overexpression impacts terminal differentiation (Takahashi et al., 1994).

Acidic and basic keratin species typically form binding partnerships within cells (Romashin et al., 2024). K16’s polymerization dimer is K6a and K6b, and loss of these two proteins results in impaired wound healing (Wong and Coulombe, 2003). Both Krt6a and Krt6b were over-expressed in the limbal compared to corneal, transitional, and conjunctival regions (Figure S10), fitting well with this cooperative. Given the many functions ascribed to Krt16, it is not surprising that its transcription is stringently controlled. As such, the transcription enhancer, nuclear factor erythroid 2-related factor 2 (Nrf2) (Yang et al., 2017), was highly expressed in the limbus vs. other regions (Figure S10), while the transcriptional repressor, Krüppel-like factor 4 (Klf4) (Szigety et al., 2020) was low in the limbus and high in the cornea and transitional zones (Figure S10).

Our data are, however, at odds with a recent report that demonstrated K16 staining confined to superficial epithelia of the mouse limbus and cornea (Wang et al., 2025). We suspect that these disparate results are due to reagent- and/or protocol-related differences, which can influence readouts. Curiously, after conducting proteomic sequencing of the corneal surface, the authors revealed low levels of K16 in the cornea compared to 13 other keratins, and by qPCR, transcripts for K16 were suppressed in the same location, results that do not align with their own immunofluorescence data. On the contrary, our data for K16 were validated by a series of molecular, biochemical, and functional approaches, while utilizing whole flat-mount specimens to provide a holistic spatial appreciation of its distribution across the ocular surface.

There is general consensus that markers such as K14, K15, p63, CD63, ABCG2, ABCB5, IFITM3, and GPHA2 are useful for identifying LESCs (Altshuler et al., 2021; Di Iorio et al., 2005; Farrelly et al., 2021; Kaplan et al., 2019; Ksander et al., 2014; Lin et al., 2023; Nakamura et al., 2008). However, they do not reliably delineate the precise location of LESCs, as their expression often overlaps with neighboring epithelia or non-epithelial cell populations (Li et al., 2024; Nureen et al., 2024; Sprogyte et al., 2024). According to our spatial transcriptomic profiling, none of these genes were differentially expressed by limbal epithelia. Nonetheless, other novel limbal markers emerged from our dataset, including Mccc2, Trim46, and Csnk2a2, which have been implicated in cancer (Chen et al., 2021; Feng et al., 2025; Ren et al., 2021). Apoe is another marker gene (Figure 3) (Li et al., 2024); however, its exclusivity for LESCs is yet to be confirmed. Noticeably, Apoe-null mice progressively develop severe corneal neovascularization (Bu et al., 2019) and suffer ocular surface dysfunction mainly due to spontaneous dyslipidemia (Módulo et al., 2012). Crtam was another upregulated limbal DEG; however, its specificity for LESC needs validation because it is also expressed by immune cells (Cortez et al., 2014), which populate the corneolimbal region (Downie et al., 2023). The only downregulated limbal DEG was Gapdh (Figure 1Ci). Given the limbus is a relatively quiescent site, it makes biological sense that a gene involved in metabolism is suppressed. Notably, similar changes in Gapdh expression were reported using scRNA-seq (Wolosin, 2024). Thus, researchers should be made aware that endogenous control genes that are often used for normalizing RNA expression, can vary across the ocular surface (Kulkarni et al., 2011).

Two decades ago, Zhou and colleagues (Zhou et al., 2006) employed laser capture microdissection to isolate basal limbal and basal corneal epithelia from fresh frozen mouse tissue sections, comparing the transcriptional signature of these two cell populations. Using the Affymetrix mouse 430 2.0 chip (Affymetrix, Santa Clara, CA), which contains the whole mouse genome, ∼100 DEGs were detected; some of the highest over-expressed genes included Epiregulin, Dachshund, and SRY, which interestingly have not been elaborated on since this report was published. At the time, this platform was regarded state of the art as the first technology to profile gene expression on-tissue, without dissociating cells from their residence. Notably, pooling of the procured cells was necessary from multiple individual mice to obtain sufficient RNA; thus, it represents a bulk sequencing assay without replicates, compared to our study that screened three independent samples in triplicate. Laser capture microdissection also risks harvesting stromal cells as well as epithelia from adjacent locations. This is a likely scenario especially when procuring limbal epithelia, which are surrounded by corneal and conjunctival counterparts. In our case, sections were stained with a pan-keratin antibody (K5) conjugated with barcoded probes to ensure only epithelial cells were sequenced. Additional assurance was provided by staining adjacent sections with H&E and immunostaining for K8/K12 to provide morphological and biochemical landmarks. These two keratin guidance markers discriminate the inner and outer limbal boundary allowing us to accurately perform a four-way comparison of DEG between cornea, transition zone, limbus, and bulbar conjunctiva (Figure 1). Although there are significant anatomical differences between the human and mouse corneolimbal landscape (Grieve et al., 2015), studies have shown that biomarkers of human LESCs, including ΔNP63, ABCB5, and NGF, also label murine basal limbal epithelia (Li et al., 2017). Likewise, there are similarities between mouse and human for K15 as a biomarker of LESC (Guo et al., 2020; Meyer-Blazejewska et al., 2010).

In the unperturbed limbus, T cells serve as niche support cells for quiescent LESCs by regulating their proliferation (Altshuler et al., 2021). In a similar vein, the normal human corneolimbal landscape is adorned with dendritic cells, monocyte/macrophages, basophils, and T cells (Li et al., 2022), mostly found within the stroma but can also be entrenched within the epithelium (Downie et al., 2023). We uncovered Nkiras1, also known as “NF-κB inhibitor-interacting Ras-like protein 1,” as an upregulated limbal DEG (Figure 1C). Nkiras1 is a tumor suppressor due to its ability to inhibit NF-κB, thereby affecting downstream circuits that regulate cell proliferation (Postler et al., 2023). NKIRAS1^+^ cells were mostly detected in the limbal stroma (Figure 3C) adjacent to the stromal blood vessels (Figure 5F) and occasionally within the limbal epithelium (Figure S4). The lobulated nuclear morphology is characteristic of neutrophils, most of which co-express LY6G (Figure 5). To our knowledge, NKIRAS1 expression has not been reported on the ocular surface, especially in neutrophils, which are numerous in the limbus, thus adding to the catalog of labels that can be used to identify these cells. Their presence ensures the limbus is protected from benign signals that would otherwise activate LESCs, thus rationalizing their position as first-line defender/responder cells. While most neutrophils were NKIRAS1^+^/LY6G^+^, a small proportion were NKIRAS1^−^/LY6G^+^ suggesting either a neutrophil subset (Carnevale et al., 2023; Matsui et al., 2020) or macrophages (Ruscitti et al., 2024) were labeled. During wound healing, we observed mobilization of Nkiras1^+^ neutrophils, suggesting their involvement in the acute-phase immune response within the limbus (Figure 5). Notably, loss of K16 leads to heightened neutrophil recruitment and amplification of type I interferon signaling (Cohen et al., 2024). As such, K16 may assist in maintaining neutrophils in check and its depletion promotes chronic inflammation and conjunctivalization that develops in LSCD. Our findings suggest that K16 may have a protective function by tempering excessive neutrophil responses and preserving limbal integrity. Another armamentarium neutrophils possess is their contribution to immune regulation by releasing neutrophil extracellular traps (NETs), which are filamentous web-like protein/nucleic acid composites that entrap and destroy microbes blocking their dissemination (Wang et al., 2024a). NETs also harbor an array of enzymes capable of degrading pro-inflammatory mediators (Schauer et al., 2014), further protecting this vital region.

To better understand the signaling circuits activated across the ocular surface, we conducted a cell-cell communication analysis, which uncovered CEACAM1 as the most enriched ligand-receptor network of signaling interactions, predominantly between limbus and conjunctiva (Figure 2). CEACAM1 is a multifaceted member of the immunoglobulin superfamily, which regulates various immune functions, including modulation of neutrophil activity, anti-microbial regulation, natural killer cell development, and stimulation of lymphocyte activity. CEACAM1 expression is elevated in peripheral blood neutrophils, suggesting a role in fine-tuning inflammatory responses and preventing chronic inflammation (Hirao et al., 2023; Lu et al., 2012; Skubitz, 2024).

IL-2 was another key signaling axis identified, arising from the limbus and directed toward conjunctiva and transitional zones (Figure 2). IL-2 plays an important role in maintaining immune homeostasis, and its dysregulation is associated with several inflammatory conditions (Mitra and Leonard, 2018). Notably, under desiccating stress, membrane-bound CD25 (the alpha chain of the IL-2 receptor) is decreased on ocular surface epithelia concomitant with an increase in its soluble species in tears. This process is mediated by matrix metalloproteinase-9 (MMP-9), which cleaves CD25 from the cell surface, thereby modulating IL-2 signaling pathways. MMP-9 activity is elevated in response to desiccating stress, suggesting its involvement in the adaptive response to environmental challenge (De Paiva et al., 2009). Neutrophils are an abundant source of MMP-9; hence these cells serve as major regulators of the primary immune response. This mechanism underscores essential limbal function in maintaining corneal integrity and immune surveillance, upon encountering environmental stressors that compromise ocular health (Zhao et al., 2023). Finally, the cornea’s enrichment of detoxification pathways underscores its critical role in protecting the ocular surface from environmental hazards, while the transitional zone appears to act as a regulatory landscape, balancing cell proliferation and differentiation to maintain ocular health.

Like all studies, our study is not exempt from limitations. The GeoMx spatial platform does not provide single-cell resolution, which restricts its ability to capture the transcriptomic profiles of individual cells. Rather it is a “mini” bulk screen that eliminates the ability to precisely explore interactions between cell types. Nonetheless, our data provide a solid platform from which DEGs can be further explored and exploited to determine epithelial cell identity. The absence of an exclusive marker of LESCs is a major impediment for their procurement and transplantation in patients with LSCD. However, caution should be exercised when interpreting data derived from animals (including mice) when comparing to humans. For example, there are obvious interspecies differences in anatomy, structure, and gene expression that are likely to affect location, function, and phenotype of these cells. Knowing their spatial location and their transcriptomic signature will provide valuable insights for precision identification and their isolation in future studies that incorporate spatial platforms at single-cell level. The genes uncovered during steady-state and their validation will heighten our understanding of pathophysiological mechanisms and lead to the discovery of novel therapeutics to treat patients with severe ocular surface disease.

Methods

Human and mouse tissue acquisition

C57BL/6 male and female wild-type mice (n = 81) were obtained from Australian BioResources (Moss Vale, Sydney, Australia). They were housed under pathogen-free conditions in temperature-controlled rooms and given ad libitum access to standard chow and water. All experimental procedures were approved by the University of New South Wales Animal Care and Ethics Committee, under protocol 23/27B. Unless specified, mice were euthanized by cervical dislocation at 8–10 weeks of age and whole eyes (left and right), consisting of the globe with conjunctiva up to the eyelids, procured. Human corneas from male and female donors were acquired according to Italian law 91/99 and guidelines from the Italian National Transplant Service (Centro Nazionale Trapianti, Rome, Italy) only after signed consent was obtained from the next of kin.

GeoMx DSP

Eyes (n = 3) were fixed in 4% paraformaldehyde (PFA) overnight. The superior nasal (SN) quadrant (Nureen et al., 2024) was dissected and paraffin-embedded (Nureen and Di Girolamo, 2024). Sections (5 μm) were mounted onto glass slides within a 35.3 mm × 14.1 mm area, dried at 58°C for 3 h, and then subjected to the GeoMx DSP workflow (Bruker Spatial Biology, San Jose, CA) (Merritt et al., 2020). Prior to selecting ROIs, serial sections were stained with hematoxylin or double immunostained for K8 and K12 to identify tissue boundaries (Nureen et al., 2024) (Figure 1). ROIs included corneal, transitional (containing some limbus and cornea), limbal, and bulbar conjunctival epithelia.

GeoMx data were analyzed using custom-R scripts. Quality control was performed using the open-source Bioconductor software package, GeoMxTools. The Seurat package was employed to perform DEA, using DESeq2 via FindAllMarkers to identify gene expression (see supplemental information for details). Multiple testing correction using the Benjamini-Hochberg procedure was applied and the significance threshold for determining DEGs corresponded to p-adj < 0.1. Analysis of cell-to-cell communication networks was performed with CellChat (v2). Spatial graphs of gene expression were generated using the GeoMx DSP software.

RT-qPCR

Whole eyes (n = 16) were dissected into corneal buttons, limbal rings, and bulbar conjunctiva. RNA was extracted from pooled samples from each region using the ISOLATE II RNA Mini Kit (Bioline, London, UK). A SuperScript IV First Strand Synthesis Kit (Invitrogen, Carlsbad, CA) was used for reverse transcription. RT-qPCR was performed using the QuantStudio 12K Flex Real-Time PCR System (Applied Biosystems, Foster City, CA) as described (Sprogyte et al., 2024). Detailed methodology is provided in supplemental information, and primers are listed in Table S1.

RNAscope workflow

Eyes (n = 4) were placed in 10% formalin overnight and then dissected and paraffin-embedded. Tissue sections (5 μm) were subjected to the RNAscope workflow in accordance with the manufacturer’s instructions (Advanced Cell Diagnostics, Minneapolis, MN). Detailed methodology is provided in supplemental information.

Immunofluorescence staining

Eyes (n = 25) were fixed in 2% PFA for 1 h at room temperature, dissected, and permeabilized with Triton X-100 (Sigma-Aldrich, St Louis, MO) and then blocked with 5% bovine serum albumin (BSA) containing 0.1% Triton X-100 diluted in phosphate-buffered saline (PBS). Tissues were incubated with primary antibodies (Table S3) for 48–72 h at 4°C and then reacted with secondary antibodies (Table S3) overnight at 4°C. Hoechst was added to counterstain nuclei, and tissues were preserved in ProLong Gold antifade DAPI mounting media (Thermo Fisher Scientific, Melbourne, Australia) prior to coverslipping. Images were acquired using confocal microscopes (Carl Zeiss, Jena, Germany). K16 intensity was quantified using ImageJ software. For other analyses, z stack images of the central cornea and one randomly selected region from each quadrant were used. Detailed methodology is provided in supplemental information.

Immunofluorescence staining in human tissue

Fresh human corneas (n = 6) were fixed in 4% PFA, then sequentially passed through a gradient of sucrose solutions, embedded in OCT compound (Sakura Finetek USA Inc., Torrance, CA), and frozen in nitrogen vapor. Cryosections were stained with anti-K16 primary antibody (Table S3) and imaged with an Eclipse Ti microscope (Nikon, Tokyo, Japan). Donor corneas were not orientated (superior/inferior, nasal/temporal) prior to sectioning and staining. Detailed methodology is provided in supplemental information. This was the only experiment conducted on human samples.

BrdU injection and associated immunodetection

BrdU (Sigma Aldrich) was intraperitoneally administered to 6-week-old C57BL/6 mice at 100 μg/g of body weight, twice a day (at 6-h intervals each day) for 4 days. Mice were euthanized 6 weeks after the last injection after which eyes were procured and stained with a rat anti-BrdU antibody (Table S3) as previously described (Nureen et al., 2024; Park et al., 2022).

Corneal debridement wound

Seven-week-old mice received a mild mechanical injury using an Algerbrush II (Katena Products, Inc., Denville, NJ) rotating burr to the epithelium within a 2-mm trephine-marked (KAI Medical, Solingen, Germany) central corneal region, as previously described (Park et al., 2019).

Statistical analysis

Unless specified, data are presented as mean ± SEM and were analyzed using GraphPad Prism (v.9.0.2) software. Statistical tests related to specific datasets can be found in supplemental information. To detect differences in mRNA fold-change, raw 2^−ΔΔCT^ values were log-transformed to minimize the possibility of skewing numerical values in each measurement variable. Unpaired Student’s t test (parametric) was used to compare differences in gene expression between limbus vs. cornea or conjunctiva in qPCR experiments and between 1-week vs. 8-week-old corneas. One-way ANOVA with Tukey’s test to correct for multiple comparisons was used to assess statistical significance where comparisons were made between cornea, limbus, and conjunctiva and between SN, superior temporal, inferior nasal, and IT quadrants. A paired Student’s t test (parametric) was used to compare parameters between wounded vs. unwounded corneas. Correlation plot of K16^+^ and BrdU^+^ cells was generated using the linear regression model. Statistical significance threshold was set at p < 0.05.

Resource availability

Lead contact

Requests for further information or resources associated with this study should be directed to Nick Di Girolamo (n.digirolamo@unsw.edu.au).

Materials availability

This study did not generate any new unique reagents.

Data and code availability

• •The datasets generated and analyzed during the current study are available in the NCBI GEO repository. The accession number for the data reported in this paper is Database: [GSE315587](GSE315587).
• •Data reported in this paper may be shared by the lead contact upon reasonable request
• •Any additional information required to re-analyze data in this report may be available upon reasonable request to the lead contact.

Acknowledgments

This work was supported by grants from the U.S Department of Defense (VR200025) and the Australian National Health and Medical Research Council (APP1156944) to N.D.G. The authors thank Dr. Behnaz Aghaei, Ms Zina Hammoudi, and Ms Jessica Leow from Invitro Technologies for their technical assistance with RNAscope assay and Professor Nicodemus Tedla, School of Biomedical Sciences, UNSW, for providing an aliquot of the Ly6G commercial antibody.

Author contributions

L.N. and A.S., investigation, data curation, methodology, formal analysis, visualization, and writing – original draft, review, and editing; S.D., data curation and investigation; V.B., resources, methodology, and writing – review and editing; S.F., resources, methodology, and writing – review and editing; D.P., resources; O.V., conceptualization, methodology, and writing – review and editing; N.D.G., conceptualization, supervision, resources, funding acquisition, methodology, and writing – original draft, review, and editing.

Declaration of interests

The authors declare no competing interests.