Cell cluster identification. Human fibrovascular membranes from patients with PDR were collected to profile the cell populations. The membranes were surgically removed and, after tissue digestion to obtain single-cell suspension, were submitted to scRNA-Seq. Raw data from 4 samples were individually processed for quality check and doublet removal (Supplemental Figure 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.172062DS1) and were then integrated in a single data set for downstream analyses. As reported in the Uniform Manifold Approximation and Projection (UMAP), we identified 6 different clusters (Figure 1A). Using canonical markers of established cell types, we identified that immune cells 1 included macrophages, expressing CD68 and IBA1, and immune cells 2 expressed high levels of CCL5 and CD3, consistent with T cells. The endothelial cells were identified based on the canonical endothelial marker CLDN5 and VWF, and the stromal cells were identified using the markers PDGFRβ, CSPG4, and ACTA2. The dot plot in Figure 1B shows the cell type classification, the percentage and average expression of specific cell markers. The cell clustering of each sample and the proportional contribution from the 4 samples to each cluster are shown in the UMAP in Figure 1C and in the bar plot in Figure 1D, respectively. Hypergeometric distribution analysis showed that sample 2 was enriched for macrophages and T cells, sample 3 was enriched for stromal cells 1 and 3, and sample 4 was enriched for T cells (Figure 1D; P < 0.001). Interestingly, endothelial cells were evenly distributed among the samples. To better investigate the molecular profile of the identified cells, we next reclustered each cell type separately.

scRNA-Seq analyses and cluster annotation. (A) Representative Uniform Manifold Approximation and Projection (UMAP) plot of the 6 different clusters revealed by Seurat analysis conducted in R Studio, and cluster identification. (B) Dot plot for common cell-specific markers, such as CD68, IBA1, CCL5, CD3, CLDN5, VWF, PDGFRB, CSPG4, and ACTA2. Cell type classification and the percentage and average expression of the specific cell markers are shown. (C and D) UMAP plot split by sample and the proportion of cells derived from each samples. Color of asterisk in D identifies which groups were enriched by hypergeometric enrichment test (P < 0.001).

Abnormal endothelial cells drive angiogenesis in PDR fibrovascular membranes. The endothelial cells were subclustered into 4 different groups using well-known enriched genes (10–12) (Figure 2A). Their relative distribution in each sample is shown in Figure 2B. The angiogenetic process is usually initiated by specialized endothelial cells, the tip cells, that sense the microenvironment, emit filopodia, and start the sprouting process, identified with the gene expression of COL4A2, MCAM, TP53I11, ESM1, and ANGPTN2 (Figure 2C). Cells expressing IGFBP3, AQP1, PLPP1, EFNB2, and ADAMTS1 were classified as stalk cells, the cells responsible for proliferating and elongating the neovascular sprout during angiogenesis. Mature endothelial cells, responsible for the most specialized functions, were identified with the expression of enriched genes such as SLC38A5, SLC7A5, and SLC32A. Immature endothelial cells were identified using the enriched genes RPS10, RPS26, RPS15, and RPS8 (Figure 2C). Immature endothelial cells are generally characterized by upregulated ribosomal gene expression and the absence of any specific endothelial-type gene expression, consistent with an intermediate phenotype (13).

Endothelial cell clustering, classification, and pathway enrichment analysis. (A and B) Representative UMAP plot of the 4 different cluster of endothelial cells is shown in A, while UMAP plot for each sample is shown in B. (C) Dot plot for the most common enriched genes for immature, stalk, tip, and mature endothelial cells. Cell type classification and the percentage and average expression of the specific cell markers are shown. (D) Dot plot of the pathway enrichment analysis for the upregulated genes of each cluster using gProfiler and PathFindR package. Only the most significant differentially expressed genes (log2FC > 1 and Padj < 0.01) were chosen for pathway enrichment analysis. The graph shows the number of genes modulated in each single pathway, the fold enrichment, and the statistical significance. (E) Dot plot showing the metabolic pathway gene expression. (F) Dot plot showing a panel of proliferation-related genes.

Pathway enrichment analysis for the upregulated genes was performed on the most significantly expressed genes (log2 fold change [log2FC > 1], adjusted P value [Padj < 0.01]). As expected, pathway-enriched genes were primarily related to a processes related to active angiogenesis in stalk and tip endothelial cells, such as blood vessel development, vasculature development, and angiogenesis as shown in the dot plots in Figure 2D. In addition, the Notch pathway, a well-known modulator of stalk cell physiology, was among the modulated pathways in this cluster. Pathways such as transporter activity regulation, vesicles, and transport across blood-brain barrier, are shown as enriched in mature endothelial cells, as expected. Immature cells showed enriched pathways related to ribosomal proteins, as expected, and glycolysis (Figure 2D). The latter pathway led us to evaluate the metabolic gene panel in all the endothelial cells. We detected a differential expression of glycolytic molecules in all the endothelial cell sublusters, a generally higher expression of these molecules compared with the tricarboxylic acid cycle genes, and a very similar expression level if compared with oxidative phosphorylation (OXPHOS) and nucleotide biosynthesis genes (Figure 2E). We observed that the metabolism was highly active in the immature endothelial cells, confirming the enrichment pathway analysis results. We also explored proliferation markers and found that they had a similar expression level among all the clusters, with the exception of immature cells, which showed low level of expression of proliferation-related genes (Figure 2F).

Pathway enrichment analysis for the downregulated genes (log2FC < –1, Padj < 0.01) is reported in the Supplemental Figure 2. As expected, genes related to angiogenesis and blood vessel development were downregulated in immature and mature endothelial cells. On the other hand, gap junction– and transporter activity–related genes were downregulated in the stalk and tip cells.

We further investigated the gene expression profile for the endothelial subclusters. Stalk cells showed the expression of expected enriched genes such as IGFBP3, ETS1, and ADAMTS1. In addition, we identified a subset of genes not physiologically expressed in stalk cells, such as VEGFC, normally involved in the formation of lymphatic vessels (14, 15), as shown in the violin plots in Figure 3A. Similarly, ESM1 and the system EFNB2/EPH4 are usually mainly expressed by tip cells, where they enhance VEGFA signaling pathway to promote the formation of new blood vessels. ANGPTN2 was found expressed in stalk cells as well, even though it is physiologically involved in modulating the migration and chemotactic activity of tip endothelial cells, by activating Rac1 (16).

Endothelial cell gene expression profile analysis. (A) Violin plot for IGFBP3, ETS1, ADAMS15, VEGFC, ESM1, EFNB2, ANGPT2, and PLAUR, which are among the most expressed genes, in each endothelial cell cluster. (B) Representative images of immunofluorescence experiments on PDR membranes for ESM1 (upper) and VEGFC (bottom). CD31 was used to stain the endothelial cells, and DAPI was used to counterstain the nuclei. Scale bar: 10 μm. Representative images from 3 independent experiments are shown. Arrowheads indicate the CD31+ESM1+ cells in the upper panel in B. They indicate VEGFC localization inside of the blood vessels in the bottom panel in B. Negative control images are shown in Supplemental Figure 2. (C) Dot plot showing the panel of genes expressed in endothelial cells forming the blood retina barrier (BRB). Cell type classification and the percentage and average expression of the specific cell markers are shown. (D) Feature plots showing BRB leakage/maintenance genes.

We next analyzed the gene expression profile of tip cells and detected expression of genes such as IGFBP3, usually highly expressed in stalk cells, and ETS1, which is normally absent in tip cells. In addition, tip cells did not express PLAUR, a gene that is physiologically tip enriched (Figure 3A).

Immunofluorescence on PDR membranes was performed to confirm the active angiogenetic process of these endothelial cells. As shown in Figure 3B, labeled with the endothelial canonical marker CD31, most of the endothelial cells coexpressed ESM1 in areas of active sprouting. The higher magnification shows colocalization of CD31 with ESM1 within these new sprouts in more detail. Next, we immunostained the membranes for VEGFC to explore the unusual expression of this growth factor as shown in scRNA-Seq. VEGFC was detected as a diffuse signal in endothelial cells, confirming its production (Figure 3B).

Next, we evaluated whether these endothelial cells showed conservation of functions normally ascribed to the retinal blood vessels from which they originated. We examined a panel of blood-retina barrier (BRB) genes normally expressed by mature endothelial cells. We found a subset of genes expressed at similar levels in tip, stalk, and mature cells, and we found an additional set of genes that were not expressed at all, as shown in the dot plot in Figure 3C. We further analyzed the expression of transcriptome specific to components of the BRB tight junctions and observed an alteration in the expression of these genes. In particular, we observed that CLDN1, CGN, and the angulin family member ILDR2 were not expressed in any of the endothelial cell clusters, while CLDN12, TJP2, OCLN, the angulin family member LSR, and the β-catenin target AXIN1 and AXIN2 were expressed at low levels in all the endothelial cell clusters. In contrast, TJP1 was highly expressed in all the endothelial cell clusters (Supplemental Figure 3). In addition, we found expression of the leakage marker PLVAP in all endothelial cells (Figure 3D), sustaining the hypothesis that those cells have lost their retina capillary phenotype. PLVAP specifically localizes to diaphragms of fenestrated endothelial cells such as the choriocapillaris (17, 18); it is typically expressed in vascular beds performing high filtration, secretion, or transendothelial transport. It is absent in endothelial cells where barrier properties are critical, such as the BRB. Interestingly, PLVAP is expressed in the BRB under pathological conditions, leading to barrier disruption (19). On the other hand, we detected upregulation of WNT pathway–related genes — such as TSPAN12, LRP5, and LRP6 — involved in the maintenance of BRB homeostasis (Figure 3D). Signaling through the frizzled class receptor 4/LDL receptor–related protein 5–6/tetraspanin 12 (FZD4/LRP5–6/TSPAN12) receptor complex is, in fact, required for developmental vascularization and BRB formation (20).

Immune cells in the fibrovascular membranes express proangiogenic genes. The importance of inflammation in the pathogenesis of DR is clear, and DR can be considered a chronic inflammatory disease (21, 22). In particular, macrophage-like cells are substantially increased on the retinal surface in human eyes with advanced DR (23–25). We focused on the immune cell clusters and identified 9 different clusters (Figure 4, A and B). To identify the cell type associated with each cluster, we relied on differentially expressed genes. CD68, CD163, CD14, and AIF1 as well as LST1 and LYZ allowed us to identify 4 macrophage clusters, and CLEC10A, CD1C, FLT3, and FCER1A were used to identify a DC cluster. The remaining 4 clusters were classified as lymphocytes. CD79A and IGKC were B cell markers; KLRB1, KLRF1, and KLRD1 identified NK cells; CCL5 and CD8A identified CD8 T cells; and FOXP3 and CD25 identified Tregs (Figure 4C).

Immune cell clustering, gene expression profile, and macrophage-expressing proangiogenetic molecule identification. (A and B) Representative UMAP plot of the 9 different cluster of immune cells is shown in A, while UMAP plot for each sample is shown in B. (C) Dot plot for the most common markers for macrophages, DCs, NK cells, and B and T lymphocytes. Cell type classification and the percentage and average expression of the specific cell markers are shown. (D) Pathway enrichment analysis for the upregulated genes of Macrophage A cluster. Only the most significant differentially expressed genes (log2FC > 1 and Padj < 0.01) were chosen for pathway enrichment analysis. The graph shows the number of genes modulated in each single pathway, the fold enrichment, and the statistical significance. (E) Dot plot of the proangiogenic cytokines, with the identification of proangiogenic SPP1+ macrophages. Cell type classification the percentage and average expression of the specific cell markers are shown. (F) Feature plot of SPP1+ cell distribution is shown in F.

It is interesting to note that CD8 T cells and NK cells were not evenly distributed among the 4 samples and were relatively underrepresented in samples 1 and 3, suggesting that these cells may not be critical to fibrovascular membranes.

Pathway enrichment analysis for the upregulated genes performed on the most expressed genes (log2FC > 1, Padj < 0.01) allowed us to identify a particular cluster of macrophages, Macro A, where chemokine mediated signaling pathways and cytokine activity were among the most modulated pathways (Figure 4D). Deeper analysis of the genes in this biological process showed the expression of chemokines known to have proangiogenetic properties, such as CXCL1, CXCL2, CXCL3, CXCL8, CCL2, and CCL3 (Figure 4E). This cluster (Macrophage A), which we labeled as proangiogenic macrophages, expressed the activated macrophage marker SPP1 and did not express microglia markers such as TMEM119 and P2RY12 (Figure 4, E and F). When compared with cluster A, the expression of SPP1 was relatively lower in clusters B and C, with dissimilar expression levels between these 2 clusters. Other immune cell clusters did not express the same panel of chemokines, but they participated in the angiogenetic process by producing other growth factors, such as ANXA2, VEGFA, and VEGFB. Ligand-receptor analysis confirmed the interaction between immune and endothelial cells. In particular, we found that macrophages communicate with the endothelial cells by transcribing molecules related to the Notch pathway, well-known to be essential for angiogenesis, in addition to the CXCL/CCL family chemokines and VEGF molecules (Supplemental Figure 4). It is worth noting that the inflammatory cells, whether macrophages or T cells, sustain a crosstalk with stromal cells as well, further supporting their potential interactions that promote the fibrovascular membranes (Supplemental Figure 4).

Pathway enrichment analysis of the clusters for the upregulated genes showed processes primarily related to specific functions of immune cells (Figure 5). Inflammatory response, phagocytosis, migration processes, granulocytes, and neutrophil chemotaxis were mainly enriched in the macrophage clusters. Enriched pathways related to MHCII activation and antigen-presenting processes were enriched in DC clusters. Pathway analysis in lymphocyte clusters revealed the activation of processes related to B cell receptor signal transduction in B lymphocytes; CD4 receptor activity in NK cells; and granzymes, cell activation, cell killing were, as expected, enriched in CD8 T cells, and lymphocyte activation–related pathways were enriched in Tregs. Pathway enrichment analysis for the downregulated genes (Supplemental Figure 5) confirmed the identity of each cluster but did not reveal particularly unexpected pathways.

Immune cell pathway enrichment analysis for the upregulated genes. Dot plot of the pathway enrichment analysis for the upregulated genes of each cluster using gProfiler and PathFindR package is shown. Only the most significant differentially expressed genes (log2FC > 1 and Padj < 0.01) were chosen for pathway enrichment analysis. The graph shows the number of genes modulated in each single pathway, the fold enrichment, and the statistical significance.

These results corroborate the important role of immune cells in fibrovascular membrane formation, showing the expression of a variety of macrophage-derived inflammatory cytokines and growth factors with proangiogenic functions.

Identifying the pericyte cluster within the stromal cells and defining their role in fibrogenesis. We next wondered if the newly formed blood vessels were covered by supporting cells that normally contribute to their function and homeostasis. We reclustered the stromal cells and identified 8 clusters (Figure 6, A and B). Based on differential gene expression, we found 1 cluster of vascular smooth muscle cells, expressing CNN2 and MYH11, but the rest of the cells were divided in 2 main groups. In the first group, we identified 2 clusters of pericytes, expressing a panel of canonical genes such as CSPG4, RGS5, KCNJ8, ABCC9, PDGFRB, and MCAM (Figure 6C), while the second group expressed myofibroblastic genes such as ACTA2, FN1, COL1A1, COL1A2, LUM, THBS2, AEBP1, MFAP5, and CTHRC1. Myofibroblasts are the cells responsible for the production of collagen and other extracellular matrix components, which are also the building blocks of scar tissue and fibrosis (Figure 6C). Interestingly, we identified an intermediate cluster between pericytes and myofibroblasts, which expressed genes overlapping with both cell types, suggesting that cells within this cluster were transitioning from their pericyte origin toward a myofibroblastic identity (Figure 6C).

Stromal cell clustering, gene expression profile, and cell trajectory inference analysis. (A and B) Representative UMAP plot of the 8 different cluster of stromal cells is shown in A, while UMAP plot for each sample is shown in B. (C) Dot plot for the most common markers for pure pericytes and myofibroblasts. Cell type classification and the percentage and average expression of the specific cell markers. (D) Trajectory inference analysis of the stromal cell clusters. The crucial node of the pseudotime analysis was identified in the transdifferentiating pericyte cluster, with 2 branches starting from this point — 1 of them toward myofibroblast lineage. (E) Dot plot of the pathway enrichment analysis for the upregulated genes of representative cluster of pericytes, transdifferentiating pericytes, and myofibroblasts, using gProfiler and PathFindR package. Only the most significant differentially expressed genes (log2FC > 1 and Padj < 0.01) were chosen for pathway enrichment analysis. The graph shows the number of genes modulated in each single pathway, the fold enrichment, and the statistical significance.

Next, we performed cell-inference trajectory analysis (or pseudotemporal ordering) to evaluate the dynamic process of transformation experienced by the cells and to identify the intermediate cell stages during this process. The analysis, reported in Figure 6D, showed that the crucial node 1 was in the transdifferentiating cluster, which we identified as the source of myofibroblasts. From this crucial node, 2 different branches emerged, describing 2 different possible dynamics. One of these branches connected with the pericyte clusters, which maintain canonical gene expression; the other branch connected to myofibroblastic clusters, with a shift in the gene expression profile, lending further support to the idea that a pericyte subset transformed toward myofibroblastic lineage (Figure 6D).

Pathway enrichment analysis for the upregulated genes reinforced this hypothesis by showing that pericytes were enriched in pathways related to physiological functions, such as blood vessel morphogenesis and development. In contrast, the transdifferentiating cluster showed pathways related to extracellular matrix components in addition to the canonical pericyte function pathways, suggesting activation of fibrogenic signaling pathways (Figure 6E). Enrichment pathway analysis in myofibroblasts showed activation of processes related to extracellular matrix components, collagen binding, extracellular matrix organization, and fibronectin binding. The pathway enrichment analysis for the downregulated genes further corroborated these results. As expected collagen binding and extracellular matrix genes were downregulated in the pericyte clusters, and genes related to angiogenesis and blood vessel development were downregulated in the myofibroblasts. The transdifferentiating pericyte cluster showed downregulation of genes related to vasculature development and the PDGFR signaling, a canonical pericyte pathway, confirming partial loss of pericyte identity (Supplemental Figure 6).

Notably, ligand receptor analysis (Supplemental Figure 4) showed an interesting cross-talk between myofibroblasts and endothelial cells through VEGFC signaling, in addition to the Notch pathway. This could be an interesting, previously unrecognized mechanism whereby VEGFC-expressing, pathological endothelial cells might promote the myofibroblast phenotype. In addition, myofibroblasts expressed molecules of the CCL/CXCL pathway and other members of the VEGF pathways as well as additional proangiogenic factors.

AEBP1 as a key marker of retinal pericyte transition toward myofibroblastic phenotype. To better clarify the molecular mechanism driving the pericyte transition toward myofibroblastic phenotype and to find molecules that might be critical for this process, we analyzed the differential gene expression in the stromal clusters. Among the top 10 most differentially expressed genes, we found THBS2, LUM, CTGF, COL8A1, VCAN, and FBLN5 significantly higher in the myofibroblasts when compared with pericytes. Most of these are well-known extracellular matrix proteins that are downstream of fibrogenic signaling pathways. Adipocyte enhancer-binding protein 1 (AEBP1), a profibrotic molecule and transcription factor, previously identified in a variety of tissues but not in preretinal fibrosis in human PDR, was also significantly upregulated in the myofibroblastic clusters when compared with the pericyte clusters, as shown in the volcano plot and in the violin plot in Figure 7A. We therefore hypothesized that AEBP1 might be one of the key molecules contributing to the pericyte change toward myofibroblastic phenotype in this context.

AEBP1 gene identification in the stromal cluster. (A) Volcano plot showing the top differentially expressed genes in the myofibroblast clusters and violin plot for AEBP1 expression in all the stromal cell clusters. (B) Immunofluorescence on PDR membranes for AEBP1 and the pericyte marker NG2. CD31 was used to visualize the endothelial cells, and DAPI was used to counterstain the nuclei. (C) Immnostaining for AEBP1, NG2, and αSMA. Scale bars: 10 μm. Arrowheads indicate NG2+AEBP1+ cells around the blood vessels; the pound symbols indicate AEBP1+ cells outside of the blood vessel area and the asterisk indicates NG2+ cells. Representative images from 3 independent experiments are shown. Negative control images are shown in Supplemental Figure 2.

Immunofluorescence of PDR membranes showed that AEBP1 was expressed in most of the NG2+ pericyte cells surrounding endothelial cells stained with CD31, reflecting the molecular change of these pericytes (Figure 7B). It is worth mentioning that some of those transitioning pericytes had slightly different morphology, appearing as elongated mesenchymal cells, highlighted with arrowheads. As expected, some of the NG2+ cells did not stain with AEBP1, suggesting that those pericytes retained their physiological gene expression and function (Figure 7B). Moreover, AEBP1 expression was also observed in other cells, not associated with the blood vessels and not expressing NG2 (Figure 7B). We cannot confirm the nature of these cells, which could be fully transdifferentiated toward myofibroblast lineage, or other cell types, unrelated to pericytes. We also stained the membranes with the myofibroblast marker αSMA and found that AEBP1+NG2+ pericytes were also αSMA+ and were localized outside of the blood vessel area; they were presumably detached pericytes undergoing transformation (Figure 7C).

To better characterize the molecular effects of the diabetic microenvironment on pericytes and corroborate the human fibrovascular PDR membranes findings, we cultured human retinal pericytes (HRP) in high-glucose medium for 24 and 48 hours and compared these with cells in osmotic control and normal media. Quantitative PCR (qPCR) analysis in HRP showed significant upregulation of profibrotic markers, such as COL1A1, COL1A2, LUM, THBS2, MFAP5, CTHRC1, and ACTA2, as early as 24 hours of high glucose. This upregulation was maintained for most of the genes at 48 hours of culture. We next used TGF-β stimulation as a positive control for its induction of the fibrotic phenotype (Figure 8A). AEBP1 gene expression was significantly upregulated as well, confirming the PDR membrane scRNA-Seq and immunofluorescence results (Figure 8A). We then analyzed the canonical pericyte genes at the same time points in culture. After 48 hours in high glucose, CSPG4, encoding for NG2, showed significant decline, suggesting the onset of pericyte cell identity loss (Figure 8B).

Molecular changes in Human Retinal Pericytes (HRP) cultured in high-glucose medium. (A and B) mRNA expression for profibrotic genes COL1A1, COL1A2, LUM, THBS2, AEBP1, CTHRC1, MFAP5, ACTA2, and FN1 (A) and pericyte markers CSPG4 and PDGFRB (B) upon 24-hour and 48-hour high glucose (30 mM) medium culture. TGF-β treatment (10 ng/mL) was used as a positive control for fibrosis induction and mannitol (24nM) was used as osmotic control. Summary of 3 independent experiments is shown. One-way ANOVA followed by Tukey’s multiple comparison test were used. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. NG, normal glucose (n = 6); Man, mannitol (n = 6); HG, high glucose (n = 6); TGF-β (n = 6).

Next, to determine whether blocking AEBP1 expression could protect pericytes against the fibrogenic effect of high glucose, we downregulated AEBP1 expression by siRNA in HRP cultured in normal or high glucose. After confirming siRNA-mediated downregulation of AEBP1 mRNA and protein expression (Figure 9A), we observed significantly decreased profibrotic gene expression in HRP cells exposed to high glucose when transfected with AEBP1 siRNA compared with cells transfected with nontargeting siRNA (scramble). Notably, we observed a substantial attenuation of COL1A1, LUM, THBS2, AEBP1, COL1A2, and CTHRC1, but not of MFAP5 and ACTA2 expression (Figure 9B). These results confirmed the involvement of AEBP1 in driving and modulating the pericyte transition toward myofibroblastic phenotype.

Molecular changes in human retinal pericytes (HRP) cultured in high-glucose medium after siAEBP1 transfection. (A) mRNA and protein expression for AEBP1 after siRNA transfection. Scale bar: 100 μm. Representative images from 3 independent experiments are shown. (B) mRNA expression for profibrotic genes COL1A1, COL1A2, LUM, THBS2, CTHRC1, MFAP5, ACTA2, and AEBP1 upon 24-hour culture in high glucose (30 mM) medium, with downregulation of endogenous AEBP1 levels by siRNA assay. Summary of 3 independent experiments is shown. Two-way ANOVA test and Mann-Whitney U test were used. *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001. NG, normal glucose (n = 6); HG, high glucose (n = 8).