scRNA-seq data identifies three TM cell clusters.

(A) Cells from limbal tissue dissected and sequenced at Columbia University (integrated B6 and 129/Sj datasets) are depicted in clusters on a UMAP space. (B) UMAP representation of subclustered trabecular meshwork (TM) containing cluster. Three distinct TM cell clusters are identified (TM1, TM2, and TM3). (C) Dendrogram illustrating the relationships among these subclusters. TM1 and TM2 are more similar to each other than to TM3. (D) Heatmap showing the signature genes for TM1, TM2, and TM3 cells. The scaled RNA expression indicates how many standard deviations a given gene’s expression is from the average expression across all examined cells (positive value is greater than the expression average and negative values are lower than expression average). (E) The scaled expression of select marker genes that robustly differentiate TM cell subsets from each other are shown with a dot plot. CE: corneal endothelium, CM: ciliary muscle, IS: iris stroma, Kera: corneal keratocytes, Scl: scleral.

Schematic of datasets and pipeline.

Each pool included limbal tissue from 8 eyes. The mouse strains and lab where data was generated are shown (Columbia University, John Lab or Duke, Stamer lab). Individual eyes were dissected to isolate a strip of limbal tissue enriched for TM cells. Limbal strips were then minced and pooled for further processing. Individual cells were sequenced using the 10X Genomics pipeline (see Methods) for sc-RNA-Seq or multiome (snRNA-seq and scATAC-seq). Schematic created with BioRender.com

Limbal cell cluster marker genes.

(A) The expression levels of various anterior segment cell type marker genes in the integrated B6 and 129 data are depicted on a dot plot. These marker genes are generally accepted to be specific to individual cell types. Epithelial cells (Krt14); neurons (Rho); endothelial cells (Egfl7); ciliary body and iris cells (Tyrp1); and immune cells (Cd52). Cluster 1 expressed multiple TM marker genes including Myoc, Acta2 (encodes α-SMA), Pitx2, and Tfap2b. Although some of the neurons may be limbal, it is possible that others are retinal contamination. (B-C) Heatmaps of differentially expressed genes across all limbal tissue cell clusters (B) and across the subclusters of cluster 1 (C).

Transcriptomic similarities between strain B6 and strain 129.

(A-B) The 129/Sj (129) and C57BL/6J (B6) strains exhibit significant overlap in marker gene expression. Heatmaps display the expression of the top 10 gene markers for each subcluster of cluster 1. The expression patterns are very similarity across strains. (C-D) UMAP representations of subclusters derived from cluster 1 separated by strain.

An independent B6 scRNA-seq dataset corroborates TM cell clusters.

(A) UMAPs of the integrated Columbia University and Duke University datasets. Cells from the different institutions occupy largely different but adjacent UMAP coordinates with some overlap, indicating batch effects. The difference was largest for TM3, but the marker genes were still conserved. This was likely driven by differences in tissue processing methods (Duke dataset, Methods). (B) Analysis of Duke University data alone independently validates our findings presented in Figure 1 from Columbia University, including the presence of 3 TM cell clusters. (C). A heatmap comparing the top 20 marker genes of each TM cell subcluster for the Columbia (C57BL/6J strain only) and Duke datasets. Overall, there is strong overlap of TM cell gene expression between datasets. Discrepancies include certain marker genes (eg. Crym) in the Columbia dataset that are not detected in the Duke dataset. These differences are associated with lower sequencing depth in the Duke data, and other technical factors/batch effects. (D-E) Dot plots showing similar TM subtype marker gene expression in the two datasets.

Comparison of our TM cells to published mouse datasets

(A) Diagrammatic representation of TM and scleral cell clusters across different mouse scRNA-seq datasets. Upper panels diagrammatically represent the current study. Lower panels represent the indicated studies. All diagrams based on marker gene expression in the cells that were clustered together. Despite some naming differences, matches for the 3 TM cell types that we identify in the current paper are evident in the previous studies. Although previous studies named additional TM subtypes, we did not subcluster to the same degree to ensure that subtype differences are robust. Previous studies had limited if any validation for some cell types. Based on our current characterization and immunofluorescent localization studies (Figures 2, S6, S8, & S9) some of the cells previously labeled TM beam cells are scleral (possibly fibroblasts). The marker gene expression for these previously named TM beam cells is highest in ‘Fibroblasts’ rather than TM cells in human scRNA-seq data, further supporting a scleral identity. In addition, cells previously named as uveal cells are a TM cell type. The cells are color-coded based on their molecular identities to match the current study. Current study subcluster identity numbers are shown. (B) Dot plot of TM cell subtype markers identified in the present study (left) using the dataset and cell type assignments (bottom) of Van Zyl et al. Based on these markers, our TM subtypes match their data as follows: TM1= JCT, TM2 = Beam A and Beam Y, and TM3 = Uveal). Scleral markers that we found are enriched in the Van Zyl Beam A and Y cells, which appear to be a mix of different cell types. (C) Expression of the Van Zyl Beam A and Y cell markers that are enriched in sclera in our combined B6 and 129 dataset (the same expression pattern was evident in all Columbia and Duke datasets in the current study). See Figure S6 for IF showing expression of these markers in Sclera but not TM.

Immunostaining demonstrates expression of specific marker genes in sclera but not TM (see Figure S5).

(A) CD34 is expressed in cells of the sclera/corneoscleral transition zone (arrows). No CD34 staining was found in the TM (> 60 sections from 4 eyes). (B) LY6C1 is expressed in the sclera. In addition to nonvascular cells (arrow), LY6C1 was expressed in vascular cells within the sclera (based on morphology and the vascular endothelial cell marker, VECAD, hash marks) and SC (yellow asterisks). LY6C1 is not found in the TM. The TM is encompassed within the dashed box while the yellow dashed line marks the inner wall of SC (see Figure S7). Scale bars = 50 µm.

Biased localization of TM cell subtypes within the TM.

(A-B) Summary diagrams showing the location of expression of TM1, TM2, and TM3 signature markers in the TM. For quantification of marker distributions, the TM was divided in half along both anterior-posterior (A) and inner-outer (B) axes (see Figure S7 localization of TM and SC). The area of the TM that was positive for each marker within each region was calculated as a percentage of the total area occupied by that marker in the entire TM on each section. The percentage of marker localization to each region was averaged across all analyzed sections. The average localization across all markers used for each TM cell subtype (see Figure S8 and S9) is shown on the diagrammatic representations of the TM using heatmaps. TM1 localization is biased towards the posterior and outer TM, whereas TM3 is biased towards the anterior and inner TM. No clear bias for TM2 was observed. A total of 286 sections were examined using immunofluorescence (IF) and in situ hybridization (ISH). Between 130-160 sections were examined for each TM cell subtype. (C-F) Representative sections for a subset of markers for each subtype used for quantification. See Figure S6 for explanation of indicated tissue orientations. (G) En face image of the TM in a 3D reconstruction of a tissue whole mount, clearly demonstrating the anterior vs posterior localization bias for TM1 (posterior bias) and TM3 (anterior bias) cell types. (H) Orthogonal 3D image that is cropped and oriented in the same planes as the tissue sections. The inner bias for α-SMA (Acta2, TM3) compared to an outer bias for MYOC (TM1) is evident. (I) More spatially refined analyses were subsequently conducted using extremely high-quality sections (50-60 per marker) by dividing the TM into 8 smaller zones (see Figure S7. In this more refined study, TM2 cells have a biased localization to a posterior and central region of TM (most enriched in zones 2 and 5, see Figure S7. (J-K) The TM subtype distributions were statistically compared across the combined zones (in parentheses) as indicated. ANOVA followed by Tukey’s honestly significant difference test. Ant: anterior TM, Pos: posterior TM, Inn: inner TM, Out: outer TM, SC: Schlemm’s canal. Grey dotted lines mark the TM zones, yellow dotted lines mark inner wall of Schlemm’s canal. All scale bars: 50 µm.

Assessing zonal distribution of TM cell subtypes.

(A) Violin plot showing expression of signature genes that were used to localize TM cell subtypes. (B-C) Diagram of flat-mounted anterior segment (B) and hematoxylin and eosin-stained sagittal section (C) indicating the X, Y and Z axis as used throughout this paper as previously reported (Kizhatil et al., 2014). (D-E) Schlemm’s canal (SC) inner wall was identified based on anatomy and the expression of Prox-GFP (green) and other endothelial cell markers. The pars plana (PP) and corneal endothelium (CE) were identified based on anatomy and DAPI staining. The posterior boundary of the TM begins just anterior to the PP and was typically aligned with the most posterior portion of SC. The anterior TM ends where the corneal endothelium begins. This typically aligned with the most anterior portion of SC. The outer boundary of the TM was immediately internal to the inner wall of SC. The inner-most TM was where the TM cells end abutting the anterior chamber. Using these boundaries, the anterior-posterior and inner-outer axis of the TM were measured. Based on these measurements for each eye, the TM was divided in half along both the inner-outer and anterior-posterior axes. (F) For more refined analysis, the anterior-posterior distance was then divided into thirds to generate posterior, central, and anterior regions. Each of these regions was then equally divided into their subregions as shown. Because the anterior TM is thinner, it was divided into two zones. These divisions generate 8 TM zones. All scale bars: 50 µm.

Distribution of TM cell subtypes by marker.

(A) Schematic representation of major marker distributions from Figure 2. (B) Additional markers for a given TM cell subtype had similar patterns of localization (see Figure S9 for example staining). However, these markers were not sampled deeply enough to perform statistical analysis. A: anterior TM, I: inner TM, O: outer TM, P: posterior TM.

Additional examples of TM subtype marker localization.

Although there is variability in expression between sections, the aggregate of all sections reveals the biased distributions of TM cell subtypes as summarized in Figures 2. Here we focus on TM expression. However, some of the markers are also expressed outside of the TM. For examples, Edn3 is expressed in vasculature, while both Lypd1 and Inmt are expressed at low levels in scleral and iris cells (but are significantly higher in TM cells). All scale bars: 50 µm. Markers assessed by IF: MYOC, CHIL1, CRYM, α-SMA, TFAP2B. Markers assessed by ISH: Inmt, Edn3, Lypd1. All scale bars: 50 µm.

Pathway analysis suggests differential involvement of TM cell subtypes in extracellular matrix regulation, growth factor signaling, and actin-binding

(A-B) Molecular comparison of all three TM cell subtypes to all other sequenced cells using gene ontology (GO) of differentially expressed genes (DEGs). The top five most significant pathways as well as three additional pathways of interest are shown for each indicated comparison (adjusted P values, X-axis was cutoff at P < 1E-10). GSEA analysis was also used to assess the enrichment or underrepresentation of these GO pathways with GSEA scores color coded on the GO charts above. Pathways significantly different by GSEA analysis are indicated by asterisks. Overall, TM cells have over-representation of various extracellular matrix (ECM) molecules/ pathways including collagens, glycosaminoglycans, and integrins compared to non-TM cells. Growth factor signaling is also enriched in TM cells compared to other cell types. (C-D) When comparing TM cell subtypes, TM1 cells are further enriched for ECM molecules such as collagens, glycosaminoglycans, and integrins as well as TGF-β signaling. (E-F) TM2-enriched ECM pathways include glycosaminoglycan binding and the laminin complex. TM2 cells are also enriched for receptor-ligand signaling and insulin growth factor binding. (G-H) Actin-binding and mitochondrial metabolism genes are enriched in TM3 cells. ECM = extracellular matrix. ECS = extracellular structure. ER = endoplasmic reticulum. IMM = inner mitochondrial membrane

TM cell pathway analysis.

(A-C) Additional analyses showing the most significant pathways (by gene ontology (GO) analysis. Pathways are further separated into enriched or underrepresented categories based on gene set enrichment analysis (GSEA) score. Asterisks indicate pathways that are significantly enriched or underrepresented based on GSEA score. ECM = extracellular matrix. ECS = extracellular structure.

Signaling interactions directed from TM cells to Schlemm’s canal and vascular endothelial cells differ across TM subtypes.

(A) Circos plot showing the top predicted interactions between TM cell ligands and Schlemm’s canal endothelial (SEC) target molecules. The top predicted targets in all TM cells are split by their expression across subtypes. TM1 and TM2 participate in more predicted signaling events compared to TM3. TM1-biased ligands include members of the Tgf-β and fibronectin signaling families. Various endothelial trophic factors are predominantly expressed in TM2 including endothelin 3, vascular endothelial growth factor A, and angiopoietin 1. See Table S5 for comprehensive list of ligand-target interaction data. (B-C) TM cell ligands that signal to collector channel (CC) endothelial cells (B) or blood endothelial cells (BECs, C). (D) The expression of all ligands predicted to have signaling interactions with SECs, CCs, or BECs, are displayed across TM cell subtypes (Dot plot).

TM-TM and SEC-TM signaling.

(A-C) Analyses showing the top predicted interactions between TM and Schlemm’s canal endothelial (SEC) ligands and individual TM cell target molecules (Circos plots). The three plots show target molecules in each individual TM cell subtype respectively. (D) The expression of all ligands predicted to have signaling interactions with TM cell subtypes (Dot plot).

Mouse TM subtype marker expression in human TM cells.

(A-C) Dot plots of the expression of the top 11 marker genes for each of our mouse TM cell subtypes that have human orthologues expressed in the Van Zyl et al human dataset. TM1 marker genes generally have higher expression in the human JCT cluster than in any other cluster of cells. Other than EDN3, TM2 and TM3 marker genes do not show a strong enrichment for the annotated TM cells compared to other cells. Among TM cells, TM2 marker genes are more enriched in beam A and beam B cells. TM3 markers are mostly expressed in relatively equal levels across all the human TM cells. Schwalbe’s line cells are also similar to TM. The anterior insertion of the TM is near Schwalbe’s line, where TM stem cells are reported to reside. Therefore, the Schwalbe’s line cells maybe TM stem cells.

Nuclear morphology of each TM subtype.

(A-C) TM1 nuclei adjacent to the SC appear shorter than TM2 and TM3 nuclei and lack the typical elongated, endothelial-like morphology of beam cells. This shorter morphology gives TM1 nuclei a more spherical appearance. All scale bars: 50 µm. (D) To rigorously assess whether TM1 nuclei are more spherical, we analyzed their 3D shapes in whole mounts, comparing them to TM3 nuclei using the ‘Sphericity’ tool in Imaris. The data show that TM1 cells indeed have more spherical nuclei, consistent with a JCT identity. However, some TM1 cells exhibit more elongated nuclei, especially when located further from the SC. A higher sphericity score indicates a more rounded nucleus (with a completely spherical shape scoring 1. We characterized 60 randomly selected nuclei for each TM cell subtype. Groups were compared using Student’s t-test.

GWAS gene expression in limbal cells.

(A) Expression of genes implicated in risk for IOP elevation (left) and primary open-angle glaucoma (POAG, right) in humans GWAS plotted against mouse limbal cell type. There is similar GWAS gene expression across TM cell subtypes and other limbal cell types. TM1 and TM2 have a slightly higher expression of GWAS genes compared to TM3. (B-C) Pathway analysis (Gene ontology, GO) shows enrichment of biological pathways for the POAG (B) or IOP (C) genes vs other expressed genes considering all cell types (see Table S6 for full list of GWAS gene-enriched pathways). (D) Enriched biological pathways among POAG GWAS genes within the indicated cell types. (E) Cell type enrichment of lipid and carbohydrate metabolism pathway genes that are significantly associated with POAG (Khawaja et al.). ECM organ. = extracellular matrix organization, cell-sub. = cell-substrate adhesion, Hormone stimulus = cellular response to hormone stimulus, small GTPase signal transduction = regulation of small GTPase mediated signal transduction, Rho signaling = regulation of Rho protein signal transduction, GF stimulus = regulation of cellular response to growth factor stimulus, Vegf signaling = vascular endothelial growth factor receptor signaling pathway.

Transcription factors dictating gene expression in TM cells

(A) A schematic representation of the pipeline used to identify open chromatin sites and active transcription factors (TFs). Individual cells were profiled using both single nucleus (sn) RNA and ATAC sequencing (multiome). TM cell clusters were identified using the snRNA-seq data, while significantly open chromatin regions were identified using the snATAC-seq data. Active TFs were determined based on the odds ratio of TF binding motifs within these open chromatin regions. Schematic created with BioRender (see here and here). (B) UMAP of subclusters derived from TM cell containing cluster (cluster 1) in the snRNA-seq data. (C-E) Example ATAC tracts for the promotor regions of selected marker genes for TM1 (C), TM2 (D), and TM3 (E). Each of these marker genes has greater promoter accessibility in the TM cell subtype in which its RNA expression is enriched (orange box). The aligned marker gene positions are shown. Coverage = normalized ATAC signal reads in transcription start site. CS = coverage scale. (F) Correlations between RNA expression levels (snRNA-seq dataset) of each TF and the chromatin accessibility levels (snATAC-seq dataset) of their respective predicted target binding motifs across all TM cell subtypes (see Methods). Select TFs with strong positive (red) or negative correlations (blue) are named.

Analysis of ATAC dataset

(A) Overlap of marker gene expression for each TM cell subtype between single-cell (sc) and single-nucleus (sn) RNA sequencing (RNA-seq). Other than the expected technical drop out of genes (no detected expression) in the snRNA-seq dataset, there is good agreement. (B) Localization of significantly open chromatin regions (snATAC-seq). (C) Confusion matrix of clusters identified by snRNA-seq and snATAC-seq after separately clustering each data type using unbiased dimensional reduction. There is a significant correlation between gene expression (RNA-based clustering) and chromatin accessibility (ATAC-based clustering) with adjusted Rand index of 0.20 (p < 0.001, permutation test). (D) Heatmaps comparing the open chromatin score at a gene promoter (left panel, snATAC-seq) to the RNA expression (right panel, snRNA-seq). Individual genes are represented on the y-axis and individual cells are plotted on the x-axis (only TM cells included). The consistent heatmap patterns indicate a strong overlap between promoter chromatin states and RNA expression for individual genes across cell types, validating the quality of these multiome datasets.

Differential molecular responses of TM cell subtypes to an Lmx1b mutation

(A) UMAPs of limbal cells from WT and Lmx1bV265D/+ (V265D) mutant mice (both genotypes B6 background). (B) Lmx1b expression across limbal cell types. Lmx1b expression is highest in TM3 cells. (C) Pathways that are changed in all TM cell subtypes. (D-F) Top pathways that are significantly altered when comparing each indicated TM cell subtype across genotypes. Many of these pathways were not significantly changed when comparing all TM cell subtypes as a group, indicating that the V265D-induced changes are strongest within each TM cell subtype. V265D has a pronounced effect on mitochondrial pathways in TM3 cells but much less so in other TM cells. ECM = extracellular matrix. ECS = extracellular structure.

Further pathway analyses of Lmx1bV265D vs WT.

(A) Heatmaps comparing marker genes for all TM cell subtypes across genotype. (B) The top 10 V265D upregulated and downregulated biological processes for the indicated comparisons across genotype are shown (all scRNA-seq data, Adjusted P values, axis cut of at 1E-10). ECM = extracellular matrix. ECS = extracellular structure.

Nicotinamide treatment protects from IOP elevation in Lmx1b mutants.

(A) Representative photos of eyes from mice of the indicated genotypes and treatments (UNT = untreated, NAM = nicotinamide treated, 550mg/kg/day provided in drinking water). NAM treatment substantially lessens anterior chamber deepening (ACD), a sensitive indicator of IOP elevation in mice. The WT and NAM-treated mutant eyes have very shallow anterior chambers while the untreated mutant eye has an obviously deepened chamber (arrowheads). (B) Distributions of anterior chamber depth based on previously defined scoring system (Tolman et al., 2021). Groups compared by Fisher’s exact test. (C-D) Boxplots of IOP (interquartile range and median line) in WT and mutant eyes of both NAM treated and untreated groups. NAM treatment significantly lessens IOP elevation in mutants compared to untreated mutant controls. Groups compared by Two-way ANOVA followed by Tukey’s Honestly significant difference. Geno = genotype.