Integrated single-cell multiomics uncovers foundational regulatory mechanisms of lens development and pathology

Abstract

Ocular lens development entails epithelial to fiber cell differentiation, defects in which cause congenital cataracts. We report the first single-cell multiomic atlas of lens development, leveraging snRNA-seq, snATAC-seq and CUT&RUN-seq to discover previously unreported mechanisms of cell fate determination and cataract-linked regulatory networks. A comprehensive profile of cis- and trans-regulatory interactions, including for the cataract-linked transcription factor MAF, is established across a temporal trajectory of fiber cell differentiation. Furthermore, we identify an epigenetic paradigm of cellular differentiation, defined by progressive loss of the H3K27 methylation writer Polycomb repressive complex 2 (PRC2). PRC2 localizes to heterochromatin domains across master-regulator transcription factor gene bodies, suggesting it safeguards epithelial cell fate. Moreover, we demonstrate that FGF hyper-stimulation in vivo leads to MAF network activation and the emergence of novel lens cell states. Collectively, these data depict a comprehensive portrait of lens fiber cell differentiation, while defining regulatory effectors of cell identity and cataract formation.

Keywords: FGF; Lens; MAF; Multiomics; PRC2; Single-cell.

Fig. 1.

Experimental overview. (A) Schematic outlines approach employed to capture a multimodal portrait of lens development. IHC, immunohistochemistry; HCR-FISH, hybridization chain reaction-fluorescent in-situ hybridization; TF, transcription factor. (B) Workflow summarizes major steps in generating single-nuclei libraries and data processing.

Fig. 2.

Single-nucleus profiling of the chick lens. (A) UMAP plot displays ASL1 abundance, in log-normalized gene counts, in nuclei from the whole eye of chicken embryos. (B) Lens nuclei were subset and reclustered, revealing three cell states. Cells are colored by ASL1 abundance, cluster, stage or cell cycle phase. (C) Violin plots display the RNA abundance for genes related to epithelial (on left) or fiber cell identity (on right). Adjusted P-values were calculated via a Wilcoxon Rank Sum test for each cluster relative to other cell populations. (D) Pathway enrichment summarizes biological functions associated with marker genes.

Fig. 3.

An inferred temporal map of epithelial to fiber cell differentiation. (A) UMAP displays lens nuclei colored by pseudotime. (B) Heatmap displays the expression of gene markers for epithelial, intermediate or fiber identity along the differentiation trajectory.

Fig. 4.

Cis-regulatory modules of the lens differentiation program. (A) UMAP generated from weighted nearest-neighbor analysis of RNA and ATAC profiles. (B) Row-normalized ATAC signals within differentially accessible regions are displayed across the pseudotime trajectory. (C,D) Annotated distribution of total accessible (C) and temporally dynamic (D) regions across the genome. (E-H) Genome browsers display accessibility signal across the loci of marker genes for each subcluster. The gene body is indicated by blue bars, with arrows oriented towards the direction of coding sequence. Links on bottom indicate predicted looping interactions between peak regions and transcription start sites. Link colors correspond to the strength of the predicted association between accessibility and expression. Peak regions linked to transcription are highlighted. The RNA abundance for each gene is displayed in violin plots on the right. Numbers along the bottom indicate genomic coordinates.

Fig. 5.

Global chromatin footprints imparted by the trans-effectors of epithelial versus fiber cell identity. (A) Heatmap displays expression patterns of a manually curated list of transcription factor-encoding genes regulated during lens fiber cell differentiation. (B) Left violin plots display RNA abundance for genes downregulated in fiber cells. Right violin plots display the chromVAR enrichment scores for the cognate motif within peak regions. Motifs are derived from the JASPAR database. (C) As in B but displaying genes upregulated in fiber cells. *adjusted P-value<0.05; ****adjusted P<0.00005 (Wilcoxon Rank Sum test comparing epithelial versus fiber clusters).

Fig. 6.

Derivation of the MAF regulatory network. (A) Top row of row-normalized heatmap displays RNA abundance of MAF RNA transcripts across the pseudotime trajectory. Rows below display accessibility of MAF motif-containing loci linked to nearby changes in transcription. (B) RNA abundances of MAF-linked genes. (C) Genes containing MAF-linked peaks are displayed, with the MAF motif highlighted and marked with an asterisk. Accessibility signal is plotted across the loci. RNA abundance for each gene is displayed in violin plots on right.

Fig. 7.

PRC2 dynamics constitute an epigenetic program of fiber cell differentiation. (A) Nuclei across the fiber cell differentiation trajectory are ordered by pseudotime on the x-axis. The y-axis represents log-scaled gene counts, and the black line represents the expression trend. Adjusted P-values were calculated via Moran's I test. (B) FISH-HCR is used to visualize changes in localization of EZH2 and JARID2 transcripts in E4 chick embryos. Scale bar: 50 µm. (C) Immunohistochemistry performed on E12 mouse sections show loss of PRC2 members during fiber cell differentiation. Scale bar: 50 µm. (D) H3K27me3 is lost from differentiating mouse fiber cells. Scale bars: 50 µm (bar in C applies to left image); 15 µm (right). (E) Schematic summarizing PRC2 dynamics.

Fig. 8.

Localization of PRC2 members and histone modifications in the E4.5 chicken lens. CUT&RUN-seq was performed for EZH2, JARID2, H3K27me3, H3K4me3 and IgG control. (A) Average signal was plotted for the top 1000 peaks for EZH2, JARID2 and H3K27me3, revealing high colocalization between targets. (B) The genomic distribution of loci co-bound by EZH2 and JARID2. (C) Top 25 genes displaying the highest signal for EZH2 and JARID2 are displayed. Duplicated genes are highlighted in red. (D) Bubble chart summarizes pathway enrichment results for top 200 genes bound by EZH2 and JARID2. (E) Genome browsers display RPGC-normalized (reads per genomic content-normalized) signal for CUT&RUN-seq targets and aggregate snATAC-seq signal. (F) Genome browsers display RPGC-normalized CUT&RUN-seq signal. Signal ranges are displayed in the top left of each track.

Fig. 9.

The FGF2-treated lens. (A) UMAP summarizes lens nuclei transcriptomes using chicken eyes 6 h after retinectomy±FGF2-coated beads. (B) UMAP displays nuclei from the retinectomy-only (control) sample, colored by the predicted cell state when annotated against the development reference dataset. (C) As in B, for FGF2-treated samples. (D,E) Feature plots display abundance of TFAP2A (E) and ASL1 (E); plots are split by condition. (F) The proportion of nuclei assigned to each cell state for each E4 condition. Red crosses indicate mean proportions. (G) Violin plots display the distribution of abundance for DEGs. **adjusted P-value<0.005; ****adjusted P-value<0.00005 (Wilcoxon Rank Sum test comparing retinectomy versus FGF2-treated lens nuclei). (H) Row-normalized heatmap displays DEGs, calculated via pseudobulk analysis. Each column corresponds to a single embryo. (I) As in H, displaying predicted MAF regulatory targets. All genes in heatmaps have adjusted P-value<0.05.

Fig. 10.

Cataract-associated genes. (A-D) Genes identified in the current study that are documented in the CAT-MAP database. Lists encompass marker genes from the developing chicken snRNA-seq dataset (A), genes with regulatory correlation to predicted MAF-binding sites (B), genes bound by PRC2 in the developing lens (C) and genes with regulation following FGF2 hyperstimulation (D). Genes are colored according to their enrichment patterns during lens development: red, epithelial enriched; blue, intermediate enriched; green, fiber enriched; black, not significantly enriched in a cluster during development or not expressed.