Single-cell atlas of DR

After performing quality control and preprocessing of the scRNA-seq data, a total of 23,724 cells were obtained for subsequent analysis. To visualize the cell types, we applied Uniform Manifold Approximation and Projection (UMAP) plot analysis to reduce the high-dimensional data. The UMAP plot revealed distinct separation of different cell types, indicating their unique transcriptional profiles (Fig. 1A). Additionally, each cell type was characterized by examining specific marker genes for individual cell clusters. These marker genes exhibited high expression levels in their respective cell types, confirming their identity and providing evidence for the accuracy of the cell type classification (Fig. 1B). By analyzing gene expression patterns, we categorized the cells into 10 distinct types, including rod cells, cone cells, amacrine cells (AC), bipolar cells (BC), Müller cells, microglia, endothelial cells (EC), horizontal cells (HC), macrophages, and pericytes. We then assessed the distribution and prevalence of these cell subtypes. The findings indicated that the rod cell subtype had the highest abundance, followed by the BC, Müller cell, and AC cell subtypes (Fig. 1C).

Fig. 1

A The UMAP plot of scRNA-seq data, colored for the 10 cell clusters. B The dotplot heatmap of the marker genes across individual cell clusters. C Stacked bar chart representing the proportional composition of various retinal cell types at different time points (0, 2, 4, and 8 weeks). Box plot showing the total counts of each retinal cell type pooled from all time points. D Top 10 GO terms in BP, CC and MF. E Top 10 KEGG pathways

To gain insights into the functional characteristics of Müller cell subtypes, we conducted a gene differential expression analysis with Seurat. This analysis aimed to investigate the differential expression of genes between Müller cells and other cell subtypes. Subsequently, we performed a gene enrichment analysis on these identified genes, detailed in Additional file 3: Table S1. In terms of BP, there is a pronounced enrichment of genes involved in the positive regulation of cell proliferation, nervous system development, and the negative regulation of neuron projection development. Additionally, processes like lipid metabolic process, cell migration, and differentiation are also notably represented. This suggests a multifaceted role of Müller cells in maintaining retinal structure and function, potentially implicating these cells in both supportive and regulatory capacities in the neural environment. The CC analysis highlights a concentration of genes associated with various structures integral to cell communication and nutrient transport, such as the plasma membrane, neuron projection, myelin sheath, and extracellular region. The presence of genes related to the Golgi apparatus and endoplasmic reticulum membrane underscores the Müller cells' involvement in protein processing and trafficking. MF results exhibit a significant overrepresentation of genes coding for proteins with binding abilities, including protein binding, oxidoreductase activity, and ion binding activities like calcium ion binding, which are critical for cell signaling and metabolic processes. The abundance of genes associated with the extracellular matrix and growth factor binding further emphasizes the role of Müller cells in extracellular matrix remodeling and interaction with growth factors, which are crucial for retinal development and repair mechanisms (Fig. 1D). Additionally, our KEGG pathway analysis identified several significantly enriched pathways, such as metabolic pathways, pathways in cancer, and fluid shear stress and atherosclerosis (Fig. 1E).

Intercellular signaling and communication network of Müller cells

To investigate the communication mechanisms of Müller cells with other cell subtypes during the progression of DR, we employed CellChat to construct the communication network between Müller cells and others, as illustrated in Additional file 1: Fig. S1. Our analysis revealed that Müller cells displayed the highest relative signal strength in both incoming and outgoing signaling patterns as depicted in Fig. 2A. This suggests that Müller cells are central to the intercellular communication networks in the DR. Specifically, we observed that pleiotrophin (PTN), prosaposin (PSAP), and midkine (MK) were the most critical signaling pathways in the information flow initiated or received by Müller cells, signifying their essential roles in the progression of DR (Fig. 2B). In the network analysis of the cellular interaction, we observed a substantial degree of complexity and specificity in cell–cell communication. The left panel (Fig. 2C), detailing the number of interactions, shows that Müller cells, microglia, and EC are hubs of cellular interaction within the retinal microenvironment, engaging in numerous contacts with other cell types such as AC, BC, and EC. The right panel, focusing on interaction strength, further elaborated on the quality of these interactions. Here, Müller cells displayed prominent signaling connections, indicating not just frequent but also strong communicative links, which may be pivotal in maintaining retinal homeostasis and responding to pathophysiological conditions. We further presented a detailed analysis of the signaling patterns and ligand–receptor interactions within the retinal microenvironment. Figure 2D elucidates a diverse suite of signaling molecules, encompassing growth factors and inflammatory mediators. Notably, Müller cells were highlighted for their extensive expression of signaling molecules, confirming their central role as previously indicated in the network diagrams. Within these cells, we identified PTN, MK, PSAP, and VEGF as the top signaling pathways. Our focus was narrowed to investigate the pathways predominantly mediated by Müller cells within the communication network. Ptn interacting with its receptor Ncl and Mdk with receptors such as Ncl and Lrp1, were prominent in Müller cell signaling output. Finally, we collectively depict the expression patterns of key signaling molecules (PTN, PSAP, MK, and VEGF) and their receptors across various cell types in the context of DR. Across all cell types, Müller cells consistently showed high expression levels of these signaling molecules, indicating their significant role in intercellular communication within the progression of DR (Fig. 2f–I).

Fig. 2

A Scatter plot depicting the incoming and outgoing interaction strengths among various cell types. B The information flow of various signaling molecules within the retinal microenvironment in the context of diabetic retinopathy. C Network analysis of the number/strength of interactions. D The heatmap of signaling pathways related to outgoing signaling patterns and incoming signaling patterns. E The dot plot of incoming/outcoming communication patterns of Müller cells. F The violin plot displaying the expression profiles of key signaling molecules PTN, PSAP, MK, and VEGF across different cell types

Temporal dynamics and cellular heterogeneity of Müller cells

Upon analyzing single-cell RNA sequencing data from a rat model, we categorized Müller cells into six distinct subgroups to delineate their dynamic phenotypic changes during the progression of diabetic retinopathy at defined intervals of 0, 2, 4, and 8 weeks, encompassing a total of 1823 Müller cells identified in the dataset (Fig. 3A, B). The clustering of gene expression patterns across the different Müller cell clusters might indicate heterogeneity in the functional states of these cells, possibly reflecting various stages within the DR environment (Fig. 3C). Pseudo-time analysis was then performed on Müller cells (Fig. 3D). We demonstrated the trajectory is bifurcated, indicating two potential differentiation paths that these cells may assume in the DR. When cross-referenced with the temporal data (0 to 8 weeks), the result showed that the Müller cells from the early weeks (0 and 2) were predominantly positioned at the beginning of the trajectory. In contrast, cells from later weeks (4 and 8) progressively occupied the terminal end of the trajectory. This pattern indicated that Müller cells undergo transcriptional changes over time in the DR condition, possibly reflecting the progression of the disease state. The distribution of Müller cell clusters along the trajectory reveals that Müller cell 1 and 0 were predominantly found at the later stages of pseudo-time. Meanwhile, Müller cell 2 and 4 appeared to be more prevalent at the early stages of the trajectory (Fig. 3D). It becomes evident that there is considerable heterogeneity in the transcriptional profiles of Müller cell subtypes during the progression of DR. The bifurcated trajectory and the chronological distribution of cells suggested that Müller cells do not represent a monolithic population but rather undergo diverse transcriptional changes over time, possibly in response to the DR milieu.

Fig. 3

A The UMAP plot of Müller cell clusters. B Stacked bar plot depicting the proportion of each Müller cell clusters’ original identity at different time points (0 week, 2 weeks, 4 weeks, and 8 weeks). C Heatmap of characteristic gene expression across Müller cell clusters. D Monocle analysis plot of Müller cells in 2D-PCA space, with a color code for pseudo-time, orig.ident and cell type. E Scatter plots illustrating the expression of Rho and Sag, across a pseudo-time trajectory in different Müller cell clusters. F Violin plots of Müller cell marker gene expression. G Gene Ontology Enrichment Analysis of the characteristic gene of Müller cell 4

We observed an unexpected high expression of the rhodopsin (Rho) within the Müller cell 4. Rhodopsin is typically expressed exclusively in rod photoreceptors, where it serves as a photopigment essential for the capture of light signals. Its expression in Müller cell 4 was indeed surprising, as these glial cells did not normally participate directly in phototransduction. Figure 3E depicts the relative expression of Rho across pseudo-time, showing a significant expression in Müller cell 4, corroborated by the UMAP plot which localizes this expression specifically to this subpopulation. Concurrently, the analysis of rod and Müller cell markers indicated that while Müller cell 4 expresses rod-specific markers like Rho and Sag, it also maintained the expression of typical Müller cell markers such as Glul and Apoe, among others (Fig. 3E–F).

This co-expression of rod and Müller cell markers suggest a phenomenon of cellular interaction or fusion, possibly through phagocytosis of damaged rod cells by Müller cell 4. Phagocytosis is a known response by Müller cells to maintain retinal integrity under stress conditions, such as diabetic retinopathy. Given the context, it is plausible to speculated that Müller cell 4 was actively engaged in the clearance of damaged rods, thus acquiring their markers. The high expression of Rho in Müller cell 4 suggested an adaptation or response mechanism that may be unique to the DR environment. The retention of rod-specific markers within the Müller cell 4 might reflect a protective mechanism aimed at preserving retinal function in the face of diabetes-induced photoreceptor damage (Fig. 3G). Gene enrichment analysis for the characteristic genes of the Müller cell 4 provided compelling evidence that supported the observations regarding the unique role of this cell type in the progression of DR. The BP enrichment analysis revealed a significant overrepresentation of genes associated with visual perception, response to stimulus, photoreceptor cell maintenance, and phototransduction. The overrepresentation of these processes in Müller cell 4 supported the hypothesis that these cells may be incorporating material from damaged rod photoreceptors through phagocytosis.

Elucidating the role of Müller cells in the progression of DR

To elucidate the complex underpinnings of DR, meticulously dissecting the retinal cell atlas is imperative. Unfortunately, the integration of scRNA-seq with transcriptomic datasets for cross-species analyses remains an underexplored territory. Such integrative approaches are critical for a holistic understanding of the nuanced roles Müller cells play throughout the progression of DR. Addressing this gap, we undertook a differential gene expression analysis leveraging transcriptomic data from GSE160306. This rich dataset includes samples from individuals diagnosed with proliferative diabetic retinopathy with concurrent diabetic macular edema (PDR-DME) and those with NPDR, as depicted in Fig. 4A.

Fig. 4

A The volcano plot of genes differentially expressed between PDR-DME and NPDR. Differentially expressed genes identified with monocle2 were labeled in the figure. B Gene co-expression regulatory network of intersecting genes. C The UMAP plot of human scRNA-seq data, colored for the 9 cell clusters. D The dotplot of the marker genes across individual cell clusters. E UMAP Visualization of Müller cell clusters and RHO expression. F UMAP plots showing the distribution of gene expression for the genes PDC, RCVRN, ROM1, and SAG across Müller cell clusters. G UMAP plots visualizing the expression patterns of RHO and PDE6G genes individually and their co-expression within Müller cell clusters

We correlated the genes displaying notable temporal changes, as identified through mouse pseudo-time analysis (Additional file 4: Table S2), with their corresponding human orthologs. This was followed by cross-referencing these genes with the outcomes of the differential gene expression analysis, as documented in Additional file 5: Table S3. Among the genes that emerged from this cross-referencing, we noted an upregulation in PDR-DME samples for TYMP, IL33, HES1, PMEL, KRT18, CEBPD, IER3, MYOC, IGF2, PRSS56, and FOXD1. Conversely, RHO, PDE6G, and THNSL2 exhibited downregulation in PDR-DME samples (Fig. 4A).

To unveil potential regulatory connections, we established a gene co-expression regulatory network involving these intersecting genes (Fig. 4B). The RHO gene was observed to have a significant positive correlation with the PDE6G gene, exhibiting a correlation coefficient of 0.92. This substantial correlation underscored the critical role of the RHO within the phototransduction pathway of the retina. As rhodopsin, encoded by the RHO gene, is integral to the visual process, its tight co-expression with PDE6G, which encodes a key component of the rod photoreceptor cGMP phosphodiesterase, highlighted the synchrony of these genes in visual function. The expression of RHO, typically restricted to rod photoreceptors, and its strong association with PDE6G in the dataset underpinned the specificity of these genes in rod photoreceptor operations. This interaction was further validated using the STRING database (Additional file 2: Fig. S2).

To validate our findings, we analyzed scRNA-seq data derived from the MR and PR of three patients each, encompassing a total of 15,441 cells, with a particular focus on identifying a spectrum of cell subtypes, including Müller cell populations (Fig. 4C–D). Figure 4E elucidates the diverse Müller cell clusters identified in the human samples. The UMAP on the right of Fig. 4E specifically illustrated the expression pattern of the RHO gene across these Müller cell clusters. Notably, the RHO gene expression was not homogeneously distributed across all Müller cell clusters. Instead, it showed a marked elevation in a subset of cells within certain clusters, highlighted by the dashed circles (Fig. 4E).

The upregulation of RHO in these Müller cell clusters may be indicative of an engagement with phagocytic activity, particularly concerning the components of damaged rod photoreceptors in the diabetic retina. This upregulation corroborated the hypothesis previously formulated based on single-cell data from SD rats (Fig. 4E).

Translating the characteristic genes identified in the Müller cell 4 from SD rats to their human orthologs, we assessed their expression levels within human Müller cell clusters. The UMAP plots depicted the expression patterns of these orthologous genes, including PDC, RCVRN, ROM1, and SAG, across the identified human Müller cell clusters. Remarkably, they exhibited localized upregulation in the same subsets of Müller cells where RHO expression was prominent. The observed gene expression trends provided additional evidence that certain Müller cell clusters, upon exposure to diabetes-induced retinal stress, may adaptively express genes that were not part of their typical repertoire. The consistent expression patterns between the rats and human data sets lent credence to the hypothesis that Müller cells are involved in preserving retinal integrity in the face of photoreceptor damage, highlighting the potential conservatism of this response mechanism across species (Fig. 4F).

Finally, we assessed the co-expression patterns of RHO and PDE6G within Müller cell clusters. The UMAP plots in Fig. 4G elucidated the co-expression patterns of RHO and PDE6G within Müller cell clusters in DR. Specifically, the concurrent expression of RHO, a gene exclusively expressed in rod photoreceptors, alongside PDE6G, a gene critical to the phototransduction process, within Müller cells suggested an atypical yet critical cellular phenomenon. Following the assessment of RHO and PDE6G co-expression within Müller cell clusters, we proceeded to validate the interaction between these two proteins with the STRING database. The visual output from the STRING analysis, as depicted in Additional file 2: Fig. S2, confirmed the interaction between RHO and PDE6G. The co-localization of these two genes within Müller cells posited a potential role for these glial cells in the phagocytic processing of rod photoreceptor components during DR. The role of Müller cells in mitigating the effects of diabetic retinopathy (DR) could be crucial, offering valuable insights into cellular resilience and the maintenance of retinal function under pathological conditions.