November 2024
Volume 65, Issue 13
Open Access
Retinal Cell Biology  |   November 2024
Single-Cell Multiomics Profiling Reveals Heterogeneity of Müller Cells in the Oxygen-Induced Retinopathy Model
Author Affiliations & Notes
  • Xueming Yao
    Department of Ophthalmology, Tianjin Medical University General Hospital, Tianjin, China
    School of Medicine, Nankai University, Tianjin, China
  • Ziqi Li
    Department of Ophthalmology, Tianjin Medical University General Hospital, Tianjin, China
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
  • Yi Lei
    Department of Pharmacology, Tianjin Key Laboratory of Inflammation Biology, School of Basic Medical Sciences, Tianjin Medical University, Tianjin, China
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
  • Qiangyun Liu
    Department of Ophthalmology, Tianjin Medical University General Hospital, Tianjin, China
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
  • Siyue Chen
    Department of Ophthalmology, Tianjin Medical University General Hospital, Tianjin, China
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
  • Haokun Zhang
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
  • Xue Dong
    Department of Pharmacology, Tianjin Key Laboratory of Inflammation Biology, School of Basic Medical Sciences, Tianjin Medical University, Tianjin, China
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
  • Kai He
    Department of Ophthalmology, Tianjin Medical University General Hospital, Tianjin, China
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
  • Ju Guo
    Department of Ophthalmology, Tianjin Medical University General Hospital, Tianjin, China
  • Mulin Jun Li
    Department of Bioinformatics, The Province and Ministry Co-Sponsored Collaborative Innovation Center for Medical Epigenetics, School of Basic Medical Sciences, Tianjin Medical University, Tianjin, China
  • Xiaohong Wang
    Department of Pharmacology, Tianjin Key Laboratory of Inflammation Biology, School of Basic Medical Sciences, Tianjin Medical University, Tianjin, China
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
  • Hua Yan
    Department of Ophthalmology, Tianjin Medical University General Hospital, Tianjin, China
    Laboratory of Molecular Ophthalmology, Tianjin Medical University, Tianjin, China
    School of Medicine, Nankai University, Tianjin, China
  • Correspondence: Mulin Jun Li, Department of Bioinformatics, School of Basic Medical Sciences, Tianjin Medical University, No. 22, Qixiang Road, Heping District, Tianjin 300070, China; [email protected]
  • Xiaohong Wang, Department of Pharmacology, School of Basic Medical Sciences, Tianjin Medical University, No. 22, Qixiang Road, Heping District, Tianjin 300070, China; [email protected]
  • Hua Yan, Department of Ophthalmology, Tianjin Medical University General Hospital, No. 154 Anshan Road, Heping District, Tianjin 30052, China; [email protected]
  • Footnotes
     XY, ZL, and YL contributed equally to the work presented here and should therefore be regarded as equivalent authors.
Investigative Ophthalmology & Visual Science November 2024, Vol.65, 8. doi:https://doi.org/10.1167/iovs.65.13.8
  • Views
  • PDF
  • Share
  • Tools
    • Alerts
      ×
      This feature is available to authenticated users only.
      Sign In or Create an Account ×
    • Get Citation

      Xueming Yao, Ziqi Li, Yi Lei, Qiangyun Liu, Siyue Chen, Haokun Zhang, Xue Dong, Kai He, Ju Guo, Mulin Jun Li, Xiaohong Wang, Hua Yan; Single-Cell Multiomics Profiling Reveals Heterogeneity of Müller Cells in the Oxygen-Induced Retinopathy Model. Invest. Ophthalmol. Vis. Sci. 2024;65(13):8. https://doi.org/10.1167/iovs.65.13.8.

      Download citation file:


      © ARVO (1962-2015); The Authors (2016-present)

      ×
  • Supplements
Abstract

Purpose: Retinal neovascularization poses heightened risks of vision loss and blindness. Despite its clinical significance, the molecular mechanisms underlying the pathogenesis of retinal neovascularization remain elusive. This study utilized single-cell multiomics profiling in an oxygen-induced retinopathy (OIR) model to comprehensively investigate the intricate molecular landscape of retinal neovascularization.

Methods: Mice were exposed to hyperoxia to induce the OIR model, and retinas were isolated for nucleus isolation. The cellular landscape of the single-nucleus suspensions was extensively characterized through single-cell multiomics sequencing. Single-cell data were integrated with genome-wide association study (GWAS) data to identify correlations between ocular cell types and diabetic retinopathy. Cell communication analysis among cells was conducted to unravel crucial ligand–receptor signals. Trajectory analysis and dynamic characterization of Müller cells were performed, followed by integration with human retinal data for pathway analysis.

Results: The multiomics dataset revealed six major ocular cell classes, with Müller cells/astrocytes showing significant associations with proliferative diabetic retinopathy (PDR). Cell communication analysis highlighted pathways that are associated with vascular proliferation and neurodevelopment, such as Vegfa–Vegfr2, Igf1–Igf1r, Nrxn3–Nlgn1, and Efna5–Epha4. Trajectory analysis identified a subset of Müller cells expressing genes linked to photoreceptor degeneration. Multiomics data integration further unveiled positively regulated genes in OIR Müller cells/astrocytes associated with axon development and neurotransmitter transmission.

Conclusions: This study significantly advances our understanding of the intricate cellular and molecular mechanisms underlying retinal neovascularization, emphasizing the pivotal role of Müller cells. The identified pathways provide valuable insights into potential therapeutic targets for PDR, offering promising directions for further research and clinical interventions.

The prevalence of pathological neovascularization in retinal vascular diseases such as diabetic retinopathy (DR), retinopathy of prematurity (ROP), and age-related macular degeneration (AMD) underscores its significance as a shared pathological feature across various ocular conditions.1 DR, a leading cause of vision loss globally, arises from microvascular complications associated with chronic hyperglycemia, whereas ROP primarily affects premature infants subjected to fluctuating oxygen levels during critical developmental stages.2,3 Conversely, AMD predominantly afflicts the elderly population and is characterized by progressive degeneration of photoreceptors and the underlying retinal pigment epithelium in the macula.4 Despite their diverse etiologies, these retinal vascular diseases share a central theme of dysregulated vascular growth, which drives the formation of abnormal blood vessels that contribute to retinal dysfunction and vision impairment. Understanding the molecular and cellular underpinnings of pathological neovascularization across these conditions is pivotal for developing targeted therapeutic strategies. 
The oxygen-induced retinopathy (OIR) model has proven to be a valuable tool in elucidating the intricate cellular and molecular changes associated with retinal neovascularization.5,6 Combining this model with single-cell sequencing techniques has enabled researchers to unravel the heterogeneity of retinal cell populations, such as stage-specific tip cells, and to discern subtle gene expression changes during the progression of vascular pathologies.7 In other single-cell RNA sequencing (scRNA-seq) studies, subpopulations of microglia have been identified, showing high expression of necroptosis-related genes, glycolytic enzymes, and the proangiogenic factor Igf1.8,9 Although these studies have provided insights into certain characteristics of the OIR model, the existing literature remains limited in exploring the open chromatin regions and epigenetic regulators at the single-cell level within the OIR model. 
The amalgamation of various omics datasets contributes to a more thorough and nuanced comprehension of biological processes and ocular disease mechanisms. Using common scRNA-seq data as an illustration, there exists a range of approaches for data integration. Integrating multiple scRNA-seq datasets of the retina has proven to be essential in revealing rare cell types, such as microglia, and uncovering heterogeneity among bipolar cells within the macular retina. Additionally, combining scRNA-seq with single-cell ATAC sequencing (scATAC-seq) facilitates the exploration of gene regulatory networks specific to certain retinal cells.1013 Furthermore, integrating scRNA-seq with genome-wide association studies (GWAS) aids in elucidating associations among genetic variations, cell types, and gene expression related to ocular diseases.14,15 In this context, the integration of diverse datasets in pathological tissues, as opposed to normal retina, is indispensable for a comprehensive understanding of the cellular mechanisms underlying retinal neovascularization. 
This study employed a comprehensive approach referred to as single-cell multiome sequencing in an OIR model to investigate the molecular intricacies of retinal neovascularization. Through meticulous data preprocessing, quality control, and analysis, the study aimed to uncover distinct cell populations, chromatin accessibility patterns, and their implications in retinal neovascularization. Furthermore, the investigation delved into the association of identified cell types with DR using GWAS data. It also explored cell communication patterns and incorporated human single-cell multiomic data, aiming to provide a comprehensive understanding of the single-cell landscape in pathological neovascularization. 
Methods
Mice and the OIR Model
The animal experiments involving mice and the OIR model were reviewed and approved by the ethics committee of Tianjin Medical University General Hospital (TMUaMEC 2020006). To construct the mouse OIR model, postnatal day 7 (P7) pups along with their nursing mothers were exposed to a hyperoxia chamber regulated by a calibrated controller and an oxygen sensor (BioSpherix, Parish, NY), where they were exposed to 75% O2 until P12. The mice were then transferred to room air until P17. 
Nuclei Isolation
Retinas from P17 mice under the OIR model were digested using 1 mg/mL Collagenase/Dispase (10269638001; Roche, Basel, Switzerland) and 60 U/mL DNase I (4716728001; Roche) in cell culture medium at 37°C for 15 minutes. The resulting single-cell suspension was filtered through a 40-µm cell strainer and rinsed with PBS containing 2% BSA. The methods for isolating, rinsing, counting, and quality control of the nuclei from cell suspensions adhered to the 10x Genomics user guide for complex tissues (CG000375; 10x Genomics, Pleasanton, CA, USA). In short, the tissues were treated with cold NP-40 lysis buffer in the presence of protease and RNase inhibitors to release nuclei from cells. Purification of the nuclei was carried out using 0.1× lysis buffer, followed by washing with a wash buffer that included RNase inhibitor and 0.1% TWEEN 20 before the permeabilization step. After the diluted nuclei buffer was resuspended, manual hemocytometry was employed to count the nuclei to achieve a targeted nuclei recovery of 10,000 nuclei. 
Single-Nucleus RNA Sequencing and Single-Nucleus ATAC Sequencing Library Generation
The preparation of combined single-nucleus RNA sequencing (snRNA-seq) and single-nucleus ATAC sequencing (snATAC-seq) libraries was conducted following the manufacturer's guidelines using the 10x Genomics Single Cell Multiome ATAC + Gene Expression kit. Sequencing of the libraries was performed using a NovaSeq 6000 sequencing system with paired-end 150-bp reads (Illumina, San Diego, CA, USA). 
snRNA-Seq and snATAC-Seq Data Preprocessing and Quality Control
The demultiplexed fastq files for snRNA-seq and snATAC-seq were processed through the Cell Ranger ARC pipeline (version 2.0.2) from 10x Genomics, yielding barcoded count matrices for gene expression and ATAC data. The count matrices were then imported into ArchR 1.0.1,16 where barcodes shared between the snRNA-seq and snATAC-seq datasets were selected. In ArchR, we applied quality-control filters by including nuclei with RNA transcript counts ranging from 200 to 50,000, with less than 10% mitochondrial reads, transcription start site (TSS) enrichment greater than 4, and more than 1000 ATAC fragments. After quality control, we filtered simulated doublets with an automated process using the ArchR filterDoublets function. 
snRNA-Seq Data Analysis
The snRNA-seq data obtained from nuclei after quality-control filtration and automatic removal of doublets underwent analysis with Seurat 4.3.0.17 Following the normalization of gene expression counts via the NormalizeData function, scaling was executed using the ScaleData function. Subsequently, the top 20 principal components were employed for graph-based clustering at a resolution of 0.5. Identification of cluster identities was manually performed by referencing marker genes (amacrine cells: Filip1l; bipolar cells: Vsx2, Otx2, Grm6, Isl1; cone cells: Opn1sw, Arr3, Opn1mw; macrophages/microglia: Adgre1, Ptprc, C1qa; Müller cells/astrocytes: Rlbp1, Glul; rod cells: Rho, Pde6a, Nr2e3, Nrl) observed in previously published snRNA-seq studies of the human and mouse retina.14,1823 Clusters expressing canonical marker genes corresponding to distinct cell types were marked as potential doublets and eliminated. Subsequent reclustering was conducted employing identical parameters. Finally, clusters representing various subpopulations of the same cell type within the ultimate dataset were grouped for subsequent analyses. 
snATAC-Seq Data Analysis
The snATAC-seq information was assessed utilizing ArchR relying on the unique barcoded cell identities obtained from the snRNA-seq. Pseudo-bulk ATAC replicates for each cell type were generated using the addGroupCoverages function with default parameters. Following this, identification of chromatin accessibility peaks outside of blacklist regions was executed utilizing the addReproduciblePeakSet function and MACS224 for snATAC-seq peaks of individual cell types. Marker peaks were discovered employing the getMarkerFeatures function with a log2 fold change ≥ 1 and a false discovery rate ≤ 0.01, as established through Wilcoxon pairwise comparisons. Promoter peaks were defined as snATAC-seq peaks located within 2000 bp upstream or 100 bp downstream of a TSS. Peaks co-accessible to promoter peaks were revealed by the getCoAccessibility function, employing a correlation cutoff of 0.3 and resolution of 1. Predicted target genes for each snATAC-seq peak were determined using the getPeak2GeneLinks function, integrating barcode-matched RNA expression data from snRNA-seq with a correlation cutoff of 0.3 and a resolution of 1. 
Sequencing Tracks and Motif Enrichment Analysis
Chromatin accessibility sequencing data were generated utilizing the plotBrowserTrack function in ArchR and subsequently normalized based on the overall read count within the TSS regions. All datasets were aligned and annotated using the mm10 reference genome unless specified otherwise. Analysis for transcription factor (TF) motif enrichment within snATAC-seq peaks was carried out via the peakAnnoEnrichment function in ArchR, utilizing default parameters and relying on position frequency matrices sourced from the CisBP (Catalog of Inferred Sequence Binding Preferences) library. The identification of TF footprints was accomplished by utilizing the getFootprints function in ArchR. To address the Tn5 insertion bias, the footprinting signals were adjusted by subtracting the Tn5 insertion signal before visualization. 
Cell-Trait Association Analysis
To elucidate the role of our data in comprehending the pathogenesis of retinal vascular diseases, diabetic retinopathy GWAS data sourced from public databases, including FinnGen,25 were strategically obtained and transformed into the hg19 genome build with LIFTOVER.26 Following the procedures outlined in the scPagwas27 manual, the GWAS data underwent necessary data pruning using PLINK28 to save computation time and were subsequently imported into R (R Foundation for Statistical Computing, Vienna, Austria) along with cell type–annotated snRNA-seq datasets. The integration of GWAS and snRNA-seq data was achieved using the scPagwas_main function in the R package scPagwas, employing a polygenic regression model to deduce cell types related to diabetic retinopathy. To ensure the efficacy of the results, the Genes_by_pathway_kegg gene list was chosen for pathway analysis as suggested. Additionally, the biomaRt29 package facilitated the conversion of mouse genes to human genes, whereby non-matching genes were subsequently removed. 
Cell Communication Analysis
The gene expression data was imported into the R package CellChat,30 specifically selecting the CellChat database for the mouse to initiate a comprehensive analysis of cell communication. Utilizing the subsetCommunication function, we precisely extracted and evaluated the ligand–receptor associations existing between Müller cells/astrocytes and bipolar cells within our dataset and publicly available data (GSE150703) from P17 OIR model mice. Also, to bolster the credibility of our findings, we integrated our dataset with the public data and assessed the ligand–receptor associations again across retinal cells. Furthermore, we conducted an in-depth examination of communication patterns among diverse cell types, in addition to the pattern signaling genes, employing the identifyCommunicationPatterns function. 
Mice
Cre recombinase activity was induced in mice via intraperitoneal injection of 100 µL of tamoxifen (T5648; Sigma-Aldrich, St. Louis, MO, USA) dissolved in corn oil (1 mg/mL) at the specified time points. Glast-CreERT2 mice (012586; The Jackson Laboratory, Bar Harbor, ME, USA) were crossed with R26-CAG-LSL-ZsGreen mice (NM-KI-200045; Shanghai Model Organisms Center, Shanghai, China) to generate offspring for subsequent sorting experiments. 
Fluorescence Activated Cell Sorting
We utilized a previously described transgenic mouse line, Glast1-CreERT2;ZsGreenLSL/+, which expresses green fluorescent protein (GFP) specifically in retinal Müller cells. Retinas from six mice were dissociated and washed several times before being incubated in a digestion buffer, comprised of RPMI 1640 Medium (Thermo Fisher Scientific), 0.25-mg/mL DNase I, and 3-mg/mL Collagenase/Dispase, for 10 minutes at 37°C. Single cells were isolated using a 40-µm strainer, and the digestion was terminated by adding 5% BSA. Cells were washed twice and analyzed using a FACSVerse flow cytometer (BD Biosciences, Franklin Lakes, NJ, USA), and the data were processed with FlowJo 10 software. GFP+ cells were enriched by sorting with the FACSVerse flow cytometer. 
RNA Extraction and Quantitative Reverse-Transcription PCR
RNA was extracted from the harvested cells using TRIzol Reagent (101254514; Sigma-Aldrich). RNA concentration was measured using a NanoDrop ND-1000 spectrophotometer (Thermo Fisher Scientific). The extracted RNA was then reverse transcribed into complementary DNA (cDNA) using the mRNA Reverse Transcription Kit (A301, TransGen Biotech, Beijing, China) according to the manufacturer's protocol. Quantitative reverse-transcription PCR (RT-qPCR) was performed with SYBR Green Master Mix (Applied Biosystems, Waltham, MA, USA) on a real-time PCR 96-well system (QuantStudio 3; Applied Biosystems). Primer sequences used in this study are listed in Supplementary Table S1. All RT-qPCR data presented are based on three biological replicates. 
Trajectory Analysis
To explore the presence of cellular trajectories within Müller cells in the OIR model, we processed the snRNA-seq data alongside publicly available scRNA-seq data (GSE150703) using the software tool Monocle 331 to evaluate the overarching patterns of cellular differentiation. We reclustered the Müller cells/astrocytes from both OIR and control conditions based on the top 25 principal components. We implemented clustering at a suitable resolution of 0.4 with the R package clustree to prevent overclustering. Subsequently, we conducted an extensive trajectory analysis on the redefined cell types. 
Immunofluorescence Staining of Cross-Sections
The mouse eyes were fixed in 4% paraformaldehyde at room temperature for 1 hour, followed by dehydration in 30%, 20%, and 10% sucrose solutions sequentially for a total of 8 hours. Subsequently, the tissues were embedded in optimal cutting temperature (OCT) compound (Sakura, Tokyo, Japan) and snap-frozen. Serial sections, each 20 µm thick, were cut using a cryostat (CM1950; Leica Biosystems, Nussloch, Germany). The sections were brought to room temperature for 20 minutes, then washed in PBS for 10 minutes. Permeabilization was carried out in 0.3% Triton X-100 in PBS for 15 minutes. The sections were then blocked with 3% BSA and 0.3% Triton X-100 in PBS for 1 hour. Primary antibody incubation was performed at 4°C for 12 hours. Following three washes with PBS, the sections were incubated with the appropriate secondary antibodies (diluted 1:300; Jackson ImmunoResearch Laboratories, West Grove, PA, USA) for 2 hours at room temperature. Images were acquired using a confocal fluorescence microscope (LSM800; Carl Zeiss Microscopy, Oberkochen, Germany), with consistent settings used for all experiments to ensure comparability. 
Human Single-Cell Multiome Data Integration
To discern differentially expressed genes within the OIR model, we first annotated the published macular single-cell multiome data derived from healthy individuals. This dataset also encompassed snRNA-seq and snATAC-seq data and was annotated using the marker genes mentioned in the literature.18 Following this, we seamlessly integrated our mouse snRNA-seq dataset with this pre-existing data through the R package Harmony.32 Differentially expressed genes to Müller cells under two distinct conditions were subsequently identified using the FindMarkers function in Seurat. The differentially expressed genes identified within the Müller cell–specific chromatin-accessible regions of the OIR model were subjected to additional gene ontology pathway enrichment analyses utilizing the R package clusterProfiler.33 
Results
Single-Cell Multiomics Profiles of Retinas Isolated From OIR Mice
To investigate genes associated with pathological retinal neovascularization, we first established the OIR model and detected neovascularization through immunofluorescence in mice (Supplementary Fig. S1). Subsequently, we meticulously isolated mouse retinas, prepared single-nucleus suspensions, and constructed a single-cell multiomics landscape of the OIR model retina through joint snRNA-seq and snATAC-seq (Fig. 1A). Following rigorous quality control and doublet removal, a robust cell atlas comprised of a total of 6347 cells was obtained. Manual annotation of 10 clusters, identified by the R package Seurat, classified these cells into six major classes: amacrine cells, bipolar cells, cone cells, macrophages/microglia, Müller cells/astrocytes, and rod cells (Fig. 1B). Among these, rod cells constituted the majority, and macrophages/microglia had the fewest cells (Fig. 1C). Figure 1D portrays the expression levels of cell-specific marker genes in our dataset, demonstrating concordance with previous studies. Notable genes include Pde6a in rod cells, Arr3 in cone cells, and Blbp1 in Müller cells/astrocytes. Moreover, we conducted a comprehensive analysis to identify marker genes specific to the six discerned cell types within the OIR model, as delineated in Supplementary Table S2
Figure 1.
 
Overview of single-cell multiomic profiling in OIR model retinal cells. (A) Comprehensive workflow schematic from the OIR model to multiomic data analysis. It begins with the dissection of retinas, followed by preparation of single nucleus suspensions, and then continues through single-cell multiome sequencing. The subsequent analyses are depicted, including motif enrichment analysis to identify regulatory DNA motifs, peak-to-peak and peak-to-gene linkage to connect regulatory elements with their target genes, cell communication networks to illustrate interactions between cell types, trajectory inference for understanding cellular differentiation paths, and the integration of human single-cell multiome data to draw comparisons and validate findings. (B) A t-Distributed Stochastic Neighbor Embedding (t-SNE) visualization of snRNA-seq clusters. Each point represents a cell, and colors denote cell types. (C) Bar plot illustrating cell composition. Different colors represent distinct cell types, with numbers above bars indicating cell counts. (D) Dot plot showing expression levels and proportions of marker genes for major cell types. Dot colors indicate expression levels, and sizes represent the proportion of cells expressing the respective marker gene. (E) t-SNE plot utilizing snATAC-seq data for dimensionality reduction, with cell annotations derived from snRNA-seq. Each point represents a cell, and colors denote cell types. (F) Heatmap of 5820 snATAC-seq marker peaks across all retinal cell clusters. Color reflects the column Z score of normalized accessibility, with red indicating higher accessibility and blue indicating lower accessibility.
Figure 1.
 
Overview of single-cell multiomic profiling in OIR model retinal cells. (A) Comprehensive workflow schematic from the OIR model to multiomic data analysis. It begins with the dissection of retinas, followed by preparation of single nucleus suspensions, and then continues through single-cell multiome sequencing. The subsequent analyses are depicted, including motif enrichment analysis to identify regulatory DNA motifs, peak-to-peak and peak-to-gene linkage to connect regulatory elements with their target genes, cell communication networks to illustrate interactions between cell types, trajectory inference for understanding cellular differentiation paths, and the integration of human single-cell multiome data to draw comparisons and validate findings. (B) A t-Distributed Stochastic Neighbor Embedding (t-SNE) visualization of snRNA-seq clusters. Each point represents a cell, and colors denote cell types. (C) Bar plot illustrating cell composition. Different colors represent distinct cell types, with numbers above bars indicating cell counts. (D) Dot plot showing expression levels and proportions of marker genes for major cell types. Dot colors indicate expression levels, and sizes represent the proportion of cells expressing the respective marker gene. (E) t-SNE plot utilizing snATAC-seq data for dimensionality reduction, with cell annotations derived from snRNA-seq. Each point represents a cell, and colors denote cell types. (F) Heatmap of 5820 snATAC-seq marker peaks across all retinal cell clusters. Color reflects the column Z score of normalized accessibility, with red indicating higher accessibility and blue indicating lower accessibility.
Cellular gene expression disparities are governed by cis-regulatory elements, necessitating an open chromatin state for binding TFs to ensure normal functionality. Employing shared barcodes from joint multiomics, we leveraged snRNA-seq clustering results to assign cells in the snATAC-seq data into the same six groups (Fig. 1E). Subsequently, using the R package ArchR and MACS2 software, we computationally identified open chromatin regions for each cell type, totaling 62,320 peaks (Supplementary Table S3). Predominantly, these peaks were situated in intronic regions (n = 27,127), followed by distal regions (n = 16,644), with exonic regions harboring the fewest peaks (n = 4569). Then, the getMarkerFeatures function was leveraged to reveal marker peaks for each cell type (Fig. 1F, Supplementary Table S4). Some cell type–specific open chromatin regions emerged near cell type marker genes identified in the snRNA-seq data, such as Arr3 in cone cells, Isl1 in bipolar cells and C1qa in macrophages/microglia (Supplementary Fig. S2A). 
To further enrich our comprehension of the regulatory landscape and gene expression dynamics within the OIR model, we revealed co-accessible peaks with promoters (Supplementary Table S5) and integrated snRNA-seq data to establish correlations between peak accessibility and gene expression (Supplementary Fig. S2B, Supplementary Table S6). With these identified snATAC-seq peaks, motif enrichment analyses were performed to predict TFs orchestrating accessible chromatin sites in each cell type, revealing multiple cell-specific motifs (Supplementary Fig. S2C). Complementing this, we executed footprint analysis on established TFs documented in the literature—Otx2 and Gsc in rod cells, Lhx2 and Pou3f1 in Müller cells/astrocytes, and Myf5 and Ascl in bipolar cells.18,34,35 This analysis revealed that motif centers are protected from Tn5 transposition, consistent with TF occupancy (Fig. 2A), reinforcing the precision of our findings. 
Figure 2.
 
Footprinting analysis and cell trait association analysis. (A) Footprinting analysis of selected TFs across different cell types. Footprints were adjusted for Tn5 insertion bias by subtracting the Tn5 insertion signal from the footprinting signal. (B) t-SNE visualization of cell trait association results, displaying adjusted P values after background correction for trait-relevant scores of each cell. Cells with an adjusted P < 0.05 are depicted in red, and those with an adjusted P > 0.05 are shown in cyan. The PDR GWAS summary statistics used were sourced from FinnGen. (C) Dot plot illustrating the correlation between DR GWAS traits from FinnGen and various cell types identified in snRNA-seq. The dashed line represents the Bonferroni-corrected P value threshold (P = 0.05/6). Colors correspond to different cell types. (D) Ranking of trait-relevant genes based on Pearson correlation coefficients (PCC) across all individual cells. Colors represent the PCC values, with red indicating higher positive correlations and blue indicating lower or negative correlations. (E) Dot plot presenting trait-relevant pathways across six ocular cell types. Dot size indicates the log-ranked P value for each pathway, and color intensity reflects the proportion of cells within each cell type genetically influenced by a given pathway.
Figure 2.
 
Footprinting analysis and cell trait association analysis. (A) Footprinting analysis of selected TFs across different cell types. Footprints were adjusted for Tn5 insertion bias by subtracting the Tn5 insertion signal from the footprinting signal. (B) t-SNE visualization of cell trait association results, displaying adjusted P values after background correction for trait-relevant scores of each cell. Cells with an adjusted P < 0.05 are depicted in red, and those with an adjusted P > 0.05 are shown in cyan. The PDR GWAS summary statistics used were sourced from FinnGen. (C) Dot plot illustrating the correlation between DR GWAS traits from FinnGen and various cell types identified in snRNA-seq. The dashed line represents the Bonferroni-corrected P value threshold (P = 0.05/6). Colors correspond to different cell types. (D) Ranking of trait-relevant genes based on Pearson correlation coefficients (PCC) across all individual cells. Colors represent the PCC values, with red indicating higher positive correlations and blue indicating lower or negative correlations. (E) Dot plot presenting trait-relevant pathways across six ocular cell types. Dot size indicates the log-ranked P value for each pathway, and color intensity reflects the proportion of cells within each cell type genetically influenced by a given pathway.
Cell Types Associated With Diabetic Retinopathy Unveiled Through Multiomic Profiling
The OIR model is widely used in researching retinal neovascular diseases such as DR. Leveraging multiomic data from the OIR model, we aimed to identify cell types associated with DR. Initially, we obtained GWAS summary statistics for DR of Europeans population from the FinnGen database (Supplementary Table S7). Utilizing the R package scPagwas, a pathway-based polygenic regression method that integrates snRNA-seq and GWAS data to elucidate the cellular context vital for complex traits diseases, our analysis revealed significant associations between proliferative diabetic retinopathy (PDR)—a severe subtype of DR characterized by neovascularization leading to severe complications such as vitreous hemorrhage and retinal detachment—and various cell groups, including amacrine cells, bipolar cells, macrophages/microglia, and Müller cells/astrocytes (Fig. 2B). 
Further analysis of DR subtypes consistently highlighted the significance of Müller cells/astrocyte across traits (Fig. 2C). This observation was validated using DR GWAS data from another cohort (Supplementary Fig. S3, Supplementary Table S8).36 Subsequently, we prioritized genes related to PDR, revealing the top 10 genes with significant Pearson correlation coefficients, as shown in Figure 2D. These genes exhibited correlations with age-related macular degeneration (q = 0.009) and retinal atrophy (q = 0.029), according to the Enrichr37 database. Additionally, the top-ranked pathways associated with PDR were prioritized, including fructose and mannose metabolism, diabetic cardiomyopathy, and the insulin signaling pathway (Fig. 2E). 
Cell Communication Insights for Müller Cells/Astrocytes and Bipolar Cells
To investigate the pathogenesis of ocular neovascularization and identify potential therapeutic avenues, we conducted a detailed cellular communication analysis using the R package CellChat. This analysis leveraged both our snRNA-seq data and publicly available scRNA-seq data from the OIR model. This analytical approach provided a comprehensive insight into the intricate interactions among cells, placing particular emphasis on signals linked to Müller cells/astrocytes, which exhibited widespread activation in cell interactions (Fig. 3A). Subsequent scrutiny of outgoing communication patterns among diverse cell types unveiled analogous communication modes in functionally similar cells. Specifically, rod cells and cone cells, two distinct photoreceptor types, shared patterns involving signals such as NOTCH, JAM, and VTN. Similarly, bipolar and amacrine cells, integral to visual information transmission, demonstrated enrichment in a distinct pattern marked by signals such as NRXN, SLIT, and EGF (Fig. 3B). Signal pathways were computed separately using both our dataset and public data, with multiple overlapping pathways underscoring the reliability of our sequencing data (Supplementary Tables S9, S10). 
Figure 3.
 
Cell communication network in the OIR model retinal cells. (A) Quantification of significant ligand–receptor pairs between any two cell populations. The width of the edges is proportionate to the number of indicated ligand–receptor pairs. (B) Visualization of inferred outgoing communication patterns from secreting cells. This Sankey diagram illustrates the associations among inferred latent patterns, cell groups, and signaling pathways. The thickness of the flow indicates the contribution of the cell group or signaling pathway to each latent pattern. (C) Highlighted significant ligand–receptor pairs specifically between Müller cells/astrocytes and bipolar cells or among themselves. The color and size of the dots denote the calculated communication probability and P values. (D) Schematic representation of the genetic mouse model used in this study. Glast1-CreERT2 mice were crossed with ZsGreenLSL/LSL mice to generate Glast1-CreERT2;ZsGreenLSL/+ offspring, which express ZsGreen specifically in Müller cells upon tamoxifen induction. (E) Experimental timeline for the OIR model and RT-qPCR analysis. (F) RT-qPCR analysis of gene expression in retinal Müller cells isolated from room air (RA) and OIR mice. The bar graphs display the relative mRNA expression levels of Vegfa, Ncam1, Cdh2, Ngfr, and Flt1, normalized to β-actin. Results are presented as mean ± SEM (n = 3 per group). Statistical significance is indicated as follows: *P < 0.05, **P < 0.01, ***P < 0.001, compared to the RA group.
Figure 3.
 
Cell communication network in the OIR model retinal cells. (A) Quantification of significant ligand–receptor pairs between any two cell populations. The width of the edges is proportionate to the number of indicated ligand–receptor pairs. (B) Visualization of inferred outgoing communication patterns from secreting cells. This Sankey diagram illustrates the associations among inferred latent patterns, cell groups, and signaling pathways. The thickness of the flow indicates the contribution of the cell group or signaling pathway to each latent pattern. (C) Highlighted significant ligand–receptor pairs specifically between Müller cells/astrocytes and bipolar cells or among themselves. The color and size of the dots denote the calculated communication probability and P values. (D) Schematic representation of the genetic mouse model used in this study. Glast1-CreERT2 mice were crossed with ZsGreenLSL/LSL mice to generate Glast1-CreERT2;ZsGreenLSL/+ offspring, which express ZsGreen specifically in Müller cells upon tamoxifen induction. (E) Experimental timeline for the OIR model and RT-qPCR analysis. (F) RT-qPCR analysis of gene expression in retinal Müller cells isolated from room air (RA) and OIR mice. The bar graphs display the relative mRNA expression levels of Vegfa, Ncam1, Cdh2, Ngfr, and Flt1, normalized to β-actin. Results are presented as mean ± SEM (n = 3 per group). Statistical significance is indicated as follows: *P < 0.05, **P < 0.01, ***P < 0.001, compared to the RA group.
Our focus then transitioned to Müller cells/astrocytes and bipolar cells, which showed significant associations with DR in cell trait association analyses. Numerous ligand–receptor pairs linked to cell proliferation were identified through integration analysis, including Vegfa–Vegfr2, Igf1–Igf1r, and Fgf9–Fgfr1.3840 Additionally, signals implicated in neuronal development and synaptic formation, encompassing Nrxn3–Nlgn1, Ncam1–Fgfr1, and Efna5–Epha4, were observed (Fig. 3C, Supplementary Table S11).4143 Notably, Müller cells/astrocytes exclusively participate in multiple ligand–receptor pathways, actively contributing to cellular adhesion, differentiation, and migration processes through the FN1 and CDH pathways. To validate these findings, we isolated Müller cells from Glast1-CreERT2;ZsGreenLSL/+ transgenic mice subjected to either OIR or room air conditions (Figs. 3D, 3E). RT-qPCR analysis of isolated cells revealed significant upregulation of Vegfa, Ncam1, Cdh2, Nlgn1, and Fn1 in the OIR group compared to the room air control group. Additionally, several other genes exhibited an upward trend in expression, although these changes did not reach statistical significance (Fig. 3F, Supplementary Fig. S4). 
Dynamic Characterization and Temporal Ordering of Müller Cells in OIR Model
To enhance our understanding of Müller cell composition in the OIR model, we conducted a meticulous reclustering analysis of the Müller/astrocyte cell group. Recognizing their inherent heterogeneity, we categorized these cells into three distinct groups by the R package clustree (Supplementary Fig. S5): one group of astrocyte cells and two groups of Müller cells (Müller-1 and Müller-2) (Fig. 4A). Müller cells from the normal retina are more inclined to enrich in Müller-1 cells, whereas in the OIR model they are distributed across the two subtypes of Müller cells (Fig. 4B). Marker genes for each cell subtype were identified and functionally annotated using the STRING database (Fig. 4C, Supplementary Table S12).44 The results revealed that marker genes for Müller-2 cells are associated with the optic nerve, participating in neuron differentiation, organization, and development. By comparing the gene expression levels of Müller cell subtypes under different conditions, we found a significant increase in the expression of marker genes (Foxp2 and Ephb1) related to Müller-2 cells (Supplementary Fig. S6A). Additionally, the results of retinal immunofluorescence staining in Glast1-CreERT2;ZsGreenLSL/+ transgenic mice after OIR modeling confirmed the presence of these two subtypes in Müller cells (Supplementary Fig. S6). 
Figure 4.
 
Directional trajectories and positive genes of Müller cells in OIR models. (A) UMAP visualization illustrating subclusters of Müller cells/astrocytes. Each point represents a cell, and colors denote cell types. (B) UMAP visualization illustrating Müller cells/astrocytes in normal and OIR mouse retinas. Cells from normal retinas are shown in yellow, and those from OIR retinas are in blue. (C) Dot plot showing expression levels and proportions of marker genes for major cell types. Dot colors indicate expression levels, and sizes represent the proportion of cells expressing the respective marker gene. (D) Single-cell inferred pseudotime trajectory graph represented in a UMAP-based embedding. Cells are ordered along pseudotime, depicted with colors ranging from purple to yellow. (E, F) Phase portraits of genes Nrxn3 and Pde4d. The upper curve plots illustrate gene expression dynamics across various pseudotime points. The lower ridge plots display gene expression levels within Müller subtypes. (G) t-SNE plot of integrated cell groups from OIR model snRNA-seq and human snRNA-seq data. Cells are colored by their cell types. (H) Dot plot illustrating Gene Ontology (GO) enriched pathways of positive genes in the OIR model. Dot size and color correspond to gene number and –log(q value), respectively.
Figure 4.
 
Directional trajectories and positive genes of Müller cells in OIR models. (A) UMAP visualization illustrating subclusters of Müller cells/astrocytes. Each point represents a cell, and colors denote cell types. (B) UMAP visualization illustrating Müller cells/astrocytes in normal and OIR mouse retinas. Cells from normal retinas are shown in yellow, and those from OIR retinas are in blue. (C) Dot plot showing expression levels and proportions of marker genes for major cell types. Dot colors indicate expression levels, and sizes represent the proportion of cells expressing the respective marker gene. (D) Single-cell inferred pseudotime trajectory graph represented in a UMAP-based embedding. Cells are ordered along pseudotime, depicted with colors ranging from purple to yellow. (E, F) Phase portraits of genes Nrxn3 and Pde4d. The upper curve plots illustrate gene expression dynamics across various pseudotime points. The lower ridge plots display gene expression levels within Müller subtypes. (G) t-SNE plot of integrated cell groups from OIR model snRNA-seq and human snRNA-seq data. Cells are colored by their cell types. (H) Dot plot illustrating Gene Ontology (GO) enriched pathways of positive genes in the OIR model. Dot size and color correspond to gene number and –log(q value), respectively.
Subsequently, we employed trajectory analysis using Monocle 3 to elucidate the dynamic changes in these two Müller cell types (Fig. 4D). The UMAP (Uniform Manifold Approximation and Projection) plot revealed potential transition relationships between Müller-1 and Müller-2 cells. We identified driver genes displaying high likelihoods and pronounced dynamic behavior, including genes Nrxn3 and Pde4d (Figs. 4E, 4F). Both genes exhibited upregulation during the transition from Müller-1 cells to Müller-2 cells. Notably, gene Nrxn3 demonstrates significant alterations in signal during cell–cell communication analysis. Additionally, gene Pde4d has been verified to be upregulated in stress-inducible photoreceptor degeneration mouse models, as well as in samples from AMD patients.45,46 
Multiomic Integration and Pathway Analysis in Müller Cells/Astrocytes
Utilizing the R package harmony, we seamlessly integrated our multiomic data with data from 49,979 nuclei derived from normal human retinas. This integration enabled us to identify differentially expressed genes in Müller cells/astrocytes across diverse conditions (Fig. 4G). We determined that 433 genes were overexpressed in Müller cells within the OIR model compared to normal samples, with 265 genes overlapping with open chromatin regions in OIR Müller/astrocyte cells (Supplementary Table S13). Subsequently, a comprehensive gene enrichment analysis was executed on these genes, identifying 251 significantly enriched pathways (Supplementary Table S14). The majority of these pathways intricately pertain to neural signaling, involving processes such as axonogenesis and synaptic signal transmission. Representative instances include the regulation of synapse organization, neuron projection development, and axon guidance, as illustrated in Figure 4H. 
To further elucidate potential therapeutic interventions for pathological retinal neovascularization, we queried the Drug–Gene Interaction Database (DGIdb)47 using the aforementioned genes. This search identified 301 approved drugs that target 50 genes in this specific geneset (Supplementary Table S15). Among these, multiple genes involved in the ligand–receptor interactions between Müller cells/astrocyte and bipolar cells have been discerned, including Fgfr2, Flt1, Epha4, Epha5, and Fn1. Gene Flt1 encodes a member of the vascular endothelial growth factor receptor family. Specific inhibitors targeting Flt1, such as sunitinib, pazopanib, and vandetanib, show promise for prospective investigations into their efficacy in the treatment of neovascular eye diseases. 
Discussion
This study employed single-cell multiomics profiling in the OIR model to comprehensively investigate the genetic and epigenetic landscape of pathological retinal neovascularization. The integration of snRNA-seq and snATAC-seq data revealed six major cell classes and their epigenetic regulation heterogeneity, with Müller cells/astrocytes showing significant associations with DR. Differentially expressed genes linked to PDR unveiled potential candidates and pathways, including fructose and mannose metabolism and the insulin signaling pathway. Cell communication analysis elucidated intricate networks, particularly highlighting interactions involving Müller cells/astrocytes and bipolar cells, providing insights into potential therapeutic targets. Dynamic characterization of Müller cells in the OIR model uncovered distinct subtypes with marker genes indicating chronological transcriptional dynamics. Multiomic integration with normal human retinas revealed differentially expressed genes enriched in neural signaling pathways, offering potential therapeutic interventions for pathological retinal neovascularization. Overall, this study advances our understanding of the cellular and molecular complexities underlying retinal vascular proliferation. 
Müller cells, the predominant glial cells in the retina that span the entire thickness of the neural retina, extend from the inner to the outer limiting membrane.48,49 They play a crucial role in supporting neuronal function and metabolism while maintaining the blood–retina barrier by regulating retinal iron and water homeostasis.50 Additionally, Müller cells participate in the phagocytosis of fragments from retinal cells, as well as foreign substances introduced during retinal development and pathological processes.51 Notably, our bioinformatics findings from cell trait association analysis also underscored a significant association between DR and Müller cells/astrocytes. Upon scrutinizing Müller cell subtypes, we observed that Müller-2 cells, in contrast to Müller-1 cells, manifest a heightened expression of genes linked to optic nerve neuron differentiation, organization, and development. This finding suggests that Müller-2 cells potentially play a more pivotal role in the initiation of the disease. 
In the progression of DR, Müller cells experience functional impairment, activating both oxidative stress and inflammatory responses.52 Serving as the primary source of vascular endothelial growth factor (VEGF) in the retina, Müller cells contribute to vascular leakage and neovascularization due to the increased levels of VEGF induced by DR.53,54 In the results of our cell communication analysis, a comprehensive network of communication between Müller cells/astrocytes and bipolar cells revealed multiple ligand–receptor pairs associated with VEGF, such as Vegfa–Vegfr1, Vegfa–Vegfr2, and Vegfa–Vegfr1r2. The identification of such signaling pathways reinforces the potential role of Müller cells in influencing vascular changes and highlights the reliability of our data in the context of retinal neovascular diseases. 
Although anti-VEGF drugs have shown success in treating age-related macular degeneration and DR, their therapeutic efficacy is generally limited.55,56 In our analysis of cell communication results, we identified additional cell proliferation–related ligands, such as insulin-like growth factor 1 (IGF1). It has been reported that IGF1 levels are elevated in the vitreous or aqueous humor of patients with PDR, and the use of IGF1 receptor antagonists in ischemia-induced mice can reduce retinal neovascularization.57,58 However, there are currently no clinically applied drugs targeting this gene. Moreover, the integration and analysis of human and murine multiomics data draw attention to the broader impact of Müller cells beyond angiogenesis, particularly in the context of neural development and synaptic signal transmission. We then identified potential drugs targeting these pathways. In summary, our findings provide potential avenues for therapeutic interventions that go beyond conventional anti-angiogenic approaches. 
Despite the valuable insights obtained from our study, it is important to acknowledge several limitations. First, identifying a limited number of cell types through multiomic sequencing may not comprehensively represent the entire retinal cellular landscape. To address this, it is recommended to either increase the number of sequenced samples or employ cell sorting techniques to enrich for specific cell types of interest. Second, there was a challenge in effectively pinpointing cell-specific open chromatin regions. Increasing the sequencing depth in future experiments would enhance the precision of snATAC-seq data, enabling a more accurate characterization of the epigenetic landscape in the OIR model. Furthermore, our study is based on data collected at a single time point, restricting our ability to capture the precise temporal dynamics of Müller cell responses. Incorporating additional time points and more control samples in future experiments would facilitate a more nuanced exploration of how Müller cell subtypes evolve over time, providing a comprehensive understanding of their functional changes during the development of the disease. Finally, although our study identified disease-associated pathways linked to pathological retinal neovascularization, these findings lack experimental validation. It is imperative for future research efforts to prioritize experimental validation, potentially through in vitro and in vivo models, to confirm the roles of these pathways in disease pathogenesis. 
Acknowledgments
Supported by grants from the National Natural Science Foundation of China (82330031, 82122018, 82020108007) and National Key Technologies R&D Program (2021YFC2401404). 
Disclosure: X. Yao, None; Z. Li, None; Y. Lei, None; Q. Liu, None; S. Chen, None; H. Zhang, None; X. Dong, None; K. He, None; J. Guo, None; M.J. Li, None; X. Wang, None; H. Yan, None 
References
Wang Z, Liu CH, Huang S, Chen J. Wnt Signaling in vascular eye diseases. Prog Retin Eye Res. 2019; 70: 110–133. [CrossRef] [PubMed]
Mazzoli V, Zhong LH, Dang VT, Shi Y, Werstuck GH. Characterization of retinal microvascular complications and the effects of endoplasmic reticulum stress in mouse models of diabetic atherosclerosis. Invest Ophthalmol Vis Sci. 2020; 61: 49. [CrossRef] [PubMed]
Hartnett ME, Penn JS. Mechanisms and management of retinopathy of prematurity. N Engl J Med. 2012; 367: 2515–2526. [CrossRef] [PubMed]
Priya RR, Chew EY, Swaroop A. Genetic studies of age-related macular degeneration: lessons, challenges, and opportunities for disease management. Ophthalmology. 2012; 119: 2526–2536. [CrossRef] [PubMed]
Zhang SX, Ma JH, Bhatta M, Fliesler SJ, Wang JJ. The unfolded protein response in retinal vascular diseases: implications and therapeutic potential beyond protein folding. Prog Retin Eye Res. 2015; 45: 111–131. [CrossRef] [PubMed]
Zhang J, Liu Z, Wu H, et al. Irisin Attenuates Pathological Neovascularization in Oxygen-Induced Retinopathy Mice. Invest Ophthalmol Vis Sci. 2022; 63: 21. [CrossRef] [PubMed]
Zarkada G, Howard JP, Xiao X, et al. Specialized endothelial tip cells guide neuroretina vascularization and blood-retina-barrier formation. Dev Cell. 2021; 56: 2237–2251.e6. [CrossRef] [PubMed]
He C, Liu Y, Huang Z, et al. A specific RIP3+ subpopulation of microglia promotes retinopathy through a hypoxia-triggered necroptotic mechanism. Proc Natl Acad Sci USA. 2021; 118: e2023290118. [CrossRef] [PubMed]
Liu Z, Shi H, Xu J, et al. Single-cell transcriptome analyses reveal microglia types associated with proliferative retinopathy. JCI Insight. 2022; 7: e160940. [CrossRef] [PubMed]
Zhang R, Huang C, Chen Y, Li T, Pang L. Single-cell transcriptomic analysis revealing changes in retinal cell subpopulation levels and the pathways involved in diabetic retinopathy. Ann Transl Med. 2022; 10: 562. [CrossRef] [PubMed]
Chen K, Wang Y, Huang Y, et al. Cross-species scRNA-seq reveals the cellular landscape of retina and early alterations in type 2 diabetes mice. Genomics. 2023; 115: 110644. [CrossRef] [PubMed]
Finkbeiner C, Ortuño-Lizarán I, Sridhar A, Hooper M, Petter S, Reh TA. Single-cell ATAC-seq of fetal human retina and stem-cell-derived retinal organoids shows changing chromatin landscapes during cell fate acquisition. Cell Rep. 2022; 38: 110294. [CrossRef] [PubMed]
Liang Q, Cheng X, Wang J, et al. A multiomics atlas of the human retina at single-cell resolution. Cell Genom. 2023; 3: 100298. [CrossRef] [PubMed]
Menon M, Mohammadi S, Davila-Velderrain J, et al. Single-cell transcriptomic atlas of the human retina identifies cell types associated with age-related macular degeneration. Nat Commun. 2019; 10: 4902. [CrossRef] [PubMed]
Wang J, Cheng X, Liang Q, et al. Single-cell multiomics of the human retina reveals hierarchical transcription factor collaboration in mediating cell type-specific effects of genetic variants on gene regulation. Genome Biol. 2023; 24: 269. [CrossRef] [PubMed]
Granja JM, Corces MR, Pierce SE, et al. ArchR is a scalable software package for integrative single-cell chromatin accessibility analysis. Nat Genet. 2021; 53: 403–411. [CrossRef] [PubMed]
Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell. 2021; 184: 3573–3587.e9. [CrossRef] [PubMed]
Wang SK, Nair S, Li R, et al. Single-cell multiome of the human retina and deep learning nominate causal variants in complex eye diseases. Cell Genom. 2022; 2: 100164. [CrossRef] [PubMed]
Choi J, Li J, Ferdous S, Liang Q, Moffitt JR, Chen R. Spatial organization of the mouse retina at single cell resolution by MERFISH. Nat Commun. 2023; 14: 4929. [CrossRef] [PubMed]
Karademir D, Todorova V, Ebner LJA, Samardzija M, Grimm C. Single-cell RNA sequencing of the retina in a model of retinitis pigmentosa reveals early responses to degeneration in rods and cones. BMC Biol. 2022; 20: 86. [CrossRef] [PubMed]
Shekhar K, Lapan SW, Whitney IE, et al. Comprehensive classification of retinal bipolar neurons by single-cell transcriptomics. Cell. 2016; 166: 1308–1323.e30. [CrossRef] [PubMed]
Cochain C, Vafadarnejad E, Arampatzi P, et al. Single-cell RNA-seq reveals the transcriptional landscape and heterogeneity of aortic macrophages in murine atherosclerosis. Circ Res. 2018; 122: 1661–1674. [CrossRef] [PubMed]
Cowan CS, Renner M, De Gennaro M, et al. Cell types of the human retina and its organoids at single-cell resolution. Cell. 2020; 182: 1623–1640.e34. [CrossRef] [PubMed]
Zhang Y, Liu T, Meyer CA, et al. Model-based analysis of ChIP-Seq (MACS). Genome Biol. 2008; 9: R137. [CrossRef] [PubMed]
Kurki MI, Karjalainen J, Palta P, et al. FinnGen provides genetic insights from a well-phenotyped isolated population. Nature. 2023; 613: 508–518. [CrossRef] [PubMed]
Hinrichs AS, Karolchik D, Baertsch R, et al. The UCSC Genome Browser Database: update 2006. Nucleic Acids Res. 2006; 34: D590–D598. [CrossRef] [PubMed]
Ma Y, Deng C, Zhou Y, et al. Polygenic regression uncovers trait-relevant cellular contexts through pathway activation transformation of single-cell RNA sequencing data. Cell Genom. 2023; 3: 100383. [CrossRef] [PubMed]
Purcell S, Neale B, Todd-Brown K, et al. PLINK: a tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007; 81: 559–575. [CrossRef] [PubMed]
Smedley D, Haider S, Ballester B, et al. BioMart—biological queries made easy. BMC Genomics. 2009; 10: 22. [CrossRef] [PubMed]
Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell–cell communication using CellChat. Nat Commun. 2021; 12: 1088. [CrossRef] [PubMed]
Cao J, Spielmann M, Qiu X, et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature. 2019; 566: 496–502. [CrossRef] [PubMed]
Korsunsky I, Millard N, Fan J, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019; 16: 1289–1296. [CrossRef] [PubMed]
Wu T, Hu E, Xu S, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation (Cambridge). 2021; 2: 100141.
Dorgau B, Collin J, Rozanska A, et al. Deciphering the spatio-temporal transcriptional and chromatin accessibility of human retinal organoid development at the single cell level. iScience. 2024; 27: 109397. [CrossRef] [PubMed]
VandenBosch LS, Wohl SG, Wilken MS, et al. Developmental changes in the accessible chromatin, transcriptome and Ascl1-binding correlate with the loss in Müller glial regenerative potential. Sci Rep. 2020; 10: 13615. [CrossRef] [PubMed]
Pollack S, Igo RP, Jr, Jensen RA, et al. Multiethnic genome-wide association study of diabetic retinopathy using liability threshold modeling of duration of diabetes and glycemic control. Diabetes. 2019; 68: 441–456. [CrossRef] [PubMed]
Kuleshov MV, Jones MR, Rouillard AD, et al. Enrichr: a comprehensive gene set enrichment analysis web server 2016 update. Nucleic Acids Res. 2016; 44: W90–W97. [CrossRef] [PubMed]
Zeng FC, Zeng MQ, Huang L, et al. Downregulation of VEGFA inhibits proliferation, promotes apoptosis, and suppresses migration and invasion of renal clear cell carcinoma. Onco Targets Ther. 2016; 9: 2131–2141. [CrossRef] [PubMed]
Ferreira Mendes JM, de Faro Valverde L, Torres Andion Vidal M, et al. Effects of IGF-1 on proliferation, angiogenesis, tumor stem cell populations and activation of AKT and Hedgehog pathways in oral squamous cell carcinoma. Int J Mol Sci. 2020; 21: 6487. [CrossRef] [PubMed]
Gamm DM, Clark E, Capowski EE, Singh R. The role of FGF9 in the production of neural retina and RPE in a pluripotent stem cell model of early human retinal development. Am J Ophthalmol. 2019; 206: 113–131. [CrossRef] [PubMed]
Gomez AM, Traunmüller L, Scheiffele P. Neurexins: molecular codes for shaping neuronal synapses. Nat Rev Neurosci. 2021; 22: 137–151. [CrossRef] [PubMed]
Vukojevic V, Mastrandreas P, Arnold A, et al. Evolutionary conserved role of neural cell adhesion molecule-1 in memory. Transl Psychiatry. 2020; 10: 217. [CrossRef] [PubMed]
Iguchi T, Oka Y, Yasumura M, et al. Mutually repulsive EphA7-EfnA5 organize region-to-region corticopontine projection by inhibiting collateral extension. J Neurosci. 2021; 41: 4795–4808. [CrossRef] [PubMed]
Szklarczyk D, Gable AL, Nastou KC, et al. The STRING database in 2021: customizable protein-protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res. 2021; 49: D605–D612. [CrossRef] [PubMed]
Luu JC, Saadane A, Leinonen H, et al. Stress resilience-enhancing drugs preserve tissue structure and function in degenerating retina via phosphodiesterase inhibition. Proc Natl Acad Sci USA. 2023; 120: e2221045120. [CrossRef] [PubMed]
Orozco LD, Owen LA, Hofmann J, et al. A systems biology approach uncovers novel disease mechanisms in age-related macular degeneration. Cell Genom. 2023; 3: 100302. [CrossRef] [PubMed]
Freshour SL, Kiwala S, Cotto KC, et al. Integration of the Drug–Gene Interaction Database (DGIdb 4.0) with open crowdsource efforts. Nucleic Acids Res. 2021; 49: D1144–D1151. [CrossRef] [PubMed]
Chen Y, Zhang T, Zeng S, et al. Transketolase in human Müller cells is critical to resist light stress through the pentose phosphate and NRF2 pathways. Redox Biol. 2022; 54: 102379. [CrossRef] [PubMed]
Zigler JS, Jr, Sinha D. βA3/A1-crystallin: more than a lens protein. Prog Retin Eye Res. 2015; 44: 62–85. [CrossRef] [PubMed]
Reichenbach A, Bringmann A. Glia of the human retina. Glia. 2020; 68: 768–796. [CrossRef] [PubMed]
Miao Y, Zhao GL, Cheng S, Wang Z, Yang XL. Activation of retinal glial cells contributes to the degeneration of ganglion cells in experimental glaucoma. Prog Retin Eye Res. 2023; 93: 101169. [CrossRef] [PubMed]
Li X, Liu J, Hoh J, Liu J. Müller cells in pathological retinal angiogenesis. Transl Res. 2019; 207: 96–106. [CrossRef] [PubMed]
Carpi-Santos R, de Melo Reis RA, Gomes FCA, Calaza KC. Contribution of Müller cells in the diabetic retinopathy development: focus on oxidative stress and inflammation. Antioxidants (Basel). 2022; 11: 617. [CrossRef] [PubMed]
Uemura A, Fruttiger M, D'Amore PA, et al. VEGFR1 signaling in retinal angiogenesis and microinflammation. Prog Retin Eye Res. 2021; 84: 100954. [CrossRef] [PubMed]
Chatziralli I, Touhami S, Cicinelli MV, et al. Disentangling the association between retinal non-perfusion and anti-VEGF agents in diabetic retinopathy. Eye (Lond). 2022; 36: 692–703. [CrossRef] [PubMed]
Sun X, Wang X, Guo X, Wang M, Liu H. Combined use of anti-VEGF drugs before and during pars plana vitrectomy for severe proliferative diabetic retinopathy. Ophthalmol Ther. 2023; 12: 3133–3142. [CrossRef] [PubMed]
Ruberte J, Ayuso E, Navarro M, et al. Increased ocular levels of IGF-1 in transgenic mice lead to diabetes-like eye disease. J Clin Invest. 2004; 113: 1149–1157. [CrossRef] [PubMed]
Mason RH, Minaker SA, Lahaie Luna G, et al. Changes in aqueous and vitreous inflammatory cytokine levels in proliferative diabetic retinopathy: a systematic review and meta-analysis [published online ahead of print June 7, 2022]. Eye (Lond).
Figure 1.
 
Overview of single-cell multiomic profiling in OIR model retinal cells. (A) Comprehensive workflow schematic from the OIR model to multiomic data analysis. It begins with the dissection of retinas, followed by preparation of single nucleus suspensions, and then continues through single-cell multiome sequencing. The subsequent analyses are depicted, including motif enrichment analysis to identify regulatory DNA motifs, peak-to-peak and peak-to-gene linkage to connect regulatory elements with their target genes, cell communication networks to illustrate interactions between cell types, trajectory inference for understanding cellular differentiation paths, and the integration of human single-cell multiome data to draw comparisons and validate findings. (B) A t-Distributed Stochastic Neighbor Embedding (t-SNE) visualization of snRNA-seq clusters. Each point represents a cell, and colors denote cell types. (C) Bar plot illustrating cell composition. Different colors represent distinct cell types, with numbers above bars indicating cell counts. (D) Dot plot showing expression levels and proportions of marker genes for major cell types. Dot colors indicate expression levels, and sizes represent the proportion of cells expressing the respective marker gene. (E) t-SNE plot utilizing snATAC-seq data for dimensionality reduction, with cell annotations derived from snRNA-seq. Each point represents a cell, and colors denote cell types. (F) Heatmap of 5820 snATAC-seq marker peaks across all retinal cell clusters. Color reflects the column Z score of normalized accessibility, with red indicating higher accessibility and blue indicating lower accessibility.
Figure 1.
 
Overview of single-cell multiomic profiling in OIR model retinal cells. (A) Comprehensive workflow schematic from the OIR model to multiomic data analysis. It begins with the dissection of retinas, followed by preparation of single nucleus suspensions, and then continues through single-cell multiome sequencing. The subsequent analyses are depicted, including motif enrichment analysis to identify regulatory DNA motifs, peak-to-peak and peak-to-gene linkage to connect regulatory elements with their target genes, cell communication networks to illustrate interactions between cell types, trajectory inference for understanding cellular differentiation paths, and the integration of human single-cell multiome data to draw comparisons and validate findings. (B) A t-Distributed Stochastic Neighbor Embedding (t-SNE) visualization of snRNA-seq clusters. Each point represents a cell, and colors denote cell types. (C) Bar plot illustrating cell composition. Different colors represent distinct cell types, with numbers above bars indicating cell counts. (D) Dot plot showing expression levels and proportions of marker genes for major cell types. Dot colors indicate expression levels, and sizes represent the proportion of cells expressing the respective marker gene. (E) t-SNE plot utilizing snATAC-seq data for dimensionality reduction, with cell annotations derived from snRNA-seq. Each point represents a cell, and colors denote cell types. (F) Heatmap of 5820 snATAC-seq marker peaks across all retinal cell clusters. Color reflects the column Z score of normalized accessibility, with red indicating higher accessibility and blue indicating lower accessibility.
Figure 2.
 
Footprinting analysis and cell trait association analysis. (A) Footprinting analysis of selected TFs across different cell types. Footprints were adjusted for Tn5 insertion bias by subtracting the Tn5 insertion signal from the footprinting signal. (B) t-SNE visualization of cell trait association results, displaying adjusted P values after background correction for trait-relevant scores of each cell. Cells with an adjusted P < 0.05 are depicted in red, and those with an adjusted P > 0.05 are shown in cyan. The PDR GWAS summary statistics used were sourced from FinnGen. (C) Dot plot illustrating the correlation between DR GWAS traits from FinnGen and various cell types identified in snRNA-seq. The dashed line represents the Bonferroni-corrected P value threshold (P = 0.05/6). Colors correspond to different cell types. (D) Ranking of trait-relevant genes based on Pearson correlation coefficients (PCC) across all individual cells. Colors represent the PCC values, with red indicating higher positive correlations and blue indicating lower or negative correlations. (E) Dot plot presenting trait-relevant pathways across six ocular cell types. Dot size indicates the log-ranked P value for each pathway, and color intensity reflects the proportion of cells within each cell type genetically influenced by a given pathway.
Figure 2.
 
Footprinting analysis and cell trait association analysis. (A) Footprinting analysis of selected TFs across different cell types. Footprints were adjusted for Tn5 insertion bias by subtracting the Tn5 insertion signal from the footprinting signal. (B) t-SNE visualization of cell trait association results, displaying adjusted P values after background correction for trait-relevant scores of each cell. Cells with an adjusted P < 0.05 are depicted in red, and those with an adjusted P > 0.05 are shown in cyan. The PDR GWAS summary statistics used were sourced from FinnGen. (C) Dot plot illustrating the correlation between DR GWAS traits from FinnGen and various cell types identified in snRNA-seq. The dashed line represents the Bonferroni-corrected P value threshold (P = 0.05/6). Colors correspond to different cell types. (D) Ranking of trait-relevant genes based on Pearson correlation coefficients (PCC) across all individual cells. Colors represent the PCC values, with red indicating higher positive correlations and blue indicating lower or negative correlations. (E) Dot plot presenting trait-relevant pathways across six ocular cell types. Dot size indicates the log-ranked P value for each pathway, and color intensity reflects the proportion of cells within each cell type genetically influenced by a given pathway.
Figure 3.
 
Cell communication network in the OIR model retinal cells. (A) Quantification of significant ligand–receptor pairs between any two cell populations. The width of the edges is proportionate to the number of indicated ligand–receptor pairs. (B) Visualization of inferred outgoing communication patterns from secreting cells. This Sankey diagram illustrates the associations among inferred latent patterns, cell groups, and signaling pathways. The thickness of the flow indicates the contribution of the cell group or signaling pathway to each latent pattern. (C) Highlighted significant ligand–receptor pairs specifically between Müller cells/astrocytes and bipolar cells or among themselves. The color and size of the dots denote the calculated communication probability and P values. (D) Schematic representation of the genetic mouse model used in this study. Glast1-CreERT2 mice were crossed with ZsGreenLSL/LSL mice to generate Glast1-CreERT2;ZsGreenLSL/+ offspring, which express ZsGreen specifically in Müller cells upon tamoxifen induction. (E) Experimental timeline for the OIR model and RT-qPCR analysis. (F) RT-qPCR analysis of gene expression in retinal Müller cells isolated from room air (RA) and OIR mice. The bar graphs display the relative mRNA expression levels of Vegfa, Ncam1, Cdh2, Ngfr, and Flt1, normalized to β-actin. Results are presented as mean ± SEM (n = 3 per group). Statistical significance is indicated as follows: *P < 0.05, **P < 0.01, ***P < 0.001, compared to the RA group.
Figure 3.
 
Cell communication network in the OIR model retinal cells. (A) Quantification of significant ligand–receptor pairs between any two cell populations. The width of the edges is proportionate to the number of indicated ligand–receptor pairs. (B) Visualization of inferred outgoing communication patterns from secreting cells. This Sankey diagram illustrates the associations among inferred latent patterns, cell groups, and signaling pathways. The thickness of the flow indicates the contribution of the cell group or signaling pathway to each latent pattern. (C) Highlighted significant ligand–receptor pairs specifically between Müller cells/astrocytes and bipolar cells or among themselves. The color and size of the dots denote the calculated communication probability and P values. (D) Schematic representation of the genetic mouse model used in this study. Glast1-CreERT2 mice were crossed with ZsGreenLSL/LSL mice to generate Glast1-CreERT2;ZsGreenLSL/+ offspring, which express ZsGreen specifically in Müller cells upon tamoxifen induction. (E) Experimental timeline for the OIR model and RT-qPCR analysis. (F) RT-qPCR analysis of gene expression in retinal Müller cells isolated from room air (RA) and OIR mice. The bar graphs display the relative mRNA expression levels of Vegfa, Ncam1, Cdh2, Ngfr, and Flt1, normalized to β-actin. Results are presented as mean ± SEM (n = 3 per group). Statistical significance is indicated as follows: *P < 0.05, **P < 0.01, ***P < 0.001, compared to the RA group.
Figure 4.
 
Directional trajectories and positive genes of Müller cells in OIR models. (A) UMAP visualization illustrating subclusters of Müller cells/astrocytes. Each point represents a cell, and colors denote cell types. (B) UMAP visualization illustrating Müller cells/astrocytes in normal and OIR mouse retinas. Cells from normal retinas are shown in yellow, and those from OIR retinas are in blue. (C) Dot plot showing expression levels and proportions of marker genes for major cell types. Dot colors indicate expression levels, and sizes represent the proportion of cells expressing the respective marker gene. (D) Single-cell inferred pseudotime trajectory graph represented in a UMAP-based embedding. Cells are ordered along pseudotime, depicted with colors ranging from purple to yellow. (E, F) Phase portraits of genes Nrxn3 and Pde4d. The upper curve plots illustrate gene expression dynamics across various pseudotime points. The lower ridge plots display gene expression levels within Müller subtypes. (G) t-SNE plot of integrated cell groups from OIR model snRNA-seq and human snRNA-seq data. Cells are colored by their cell types. (H) Dot plot illustrating Gene Ontology (GO) enriched pathways of positive genes in the OIR model. Dot size and color correspond to gene number and –log(q value), respectively.
Figure 4.
 
Directional trajectories and positive genes of Müller cells in OIR models. (A) UMAP visualization illustrating subclusters of Müller cells/astrocytes. Each point represents a cell, and colors denote cell types. (B) UMAP visualization illustrating Müller cells/astrocytes in normal and OIR mouse retinas. Cells from normal retinas are shown in yellow, and those from OIR retinas are in blue. (C) Dot plot showing expression levels and proportions of marker genes for major cell types. Dot colors indicate expression levels, and sizes represent the proportion of cells expressing the respective marker gene. (D) Single-cell inferred pseudotime trajectory graph represented in a UMAP-based embedding. Cells are ordered along pseudotime, depicted with colors ranging from purple to yellow. (E, F) Phase portraits of genes Nrxn3 and Pde4d. The upper curve plots illustrate gene expression dynamics across various pseudotime points. The lower ridge plots display gene expression levels within Müller subtypes. (G) t-SNE plot of integrated cell groups from OIR model snRNA-seq and human snRNA-seq data. Cells are colored by their cell types. (H) Dot plot illustrating Gene Ontology (GO) enriched pathways of positive genes in the OIR model. Dot size and color correspond to gene number and –log(q value), respectively.
×
×

This PDF is available to Subscribers Only

Sign in or purchase a subscription to access this content. ×

You must be signed into an individual account to use this feature.

×