Bulk Gene Expression Deconvolution Reveals Infiltration of M2 Macrophages in Retinal Neovascularization

Abstract

Purpose: This study interrogated the transcriptional features and immune cellular landscape of the retinae of rats subjected to oxygen-induced retinopathy (OIR).

Methods: Bulk RNA sequencing was performed with retinal RNA isolated from control and OIR rats. Gene set enrichment analysis (GSEA) was undertaken to identify gene sets associated with immune responses in retinal neovascularization. Bulk gene expression deconvolution analysis by CIBERSORTx was performed to identify immune cell types involved in retinal neovascularization, followed by functional enrichment analysis of differentially expressed genes (DEGs). Protein–protein interaction analysis was performed to predict the hub genes relevant to identified immune cell types. CIBERSORTx was applied to profile immune cell types in the macula of patients with both proliferative diabetic retinopathy (PDR) and diabetic macular edema using a public RNA-seq dataset.

Results: Transcriptome analysis by GSEA revealed that the retina of OIR rats and patients with PDR is characterized by increased immunoregulatory interactions and complement cascade. Deconvolution analysis demonstrated that M2 macrophages infiltrate the retinae of OIR rats and patients with PDR. Functional enrichment analysis of DEGs in OIR rats showed that the dysregulated genes are related to leukocyte-mediated immunity and myeloid leukocyte activation. Downstream protein–protein interaction analysis revealed that several potential hub genes, including Ccl2, Itgam, and Tlr2, contribute to M2 macrophage infiltration in the ischemic retina.

Conclusions: This study highlights application of the gene expression deconvolution tool to identify immune cell types in inflammatory ocular diseases with transcriptomes, providing a new approach to assess changes in immune cell types in diseased ocular tissues.

Retinal neovascularization is a severe complication of a number of ocular diseases, most commonly in proliferative diabetic retinopathy (PDR), resulting in irreversible visual impairment or vision loss.¹ PDR has classically been known as a microvascular complication of diabetes and is characterized by vascular abnormalities.² However, recent discoveries have indicated involvement of the immune system, particularly the complement system, in maintaining retinal integrity, disease pathogenesis, and severity.^3,4 In line with these findings, a number of inflammatory cytokines and chemokines have been reported to be elevated in ocular tissues of both patients with PDR and animal models.⁵ 

The mechanism behind these immune response-based manifestations in patients with retinal neovascularization remains unclear. This may, in part, be due to the complexity of the immune responses and the various immune cell types, cytokines and chemokines, involved. In addition, it is still difficult to understand whether the immune system is activated locally in the retina or systemically, further complicating the question. However, a growing body of evidence based on the use of different experimental animal models has concluded that there is a common feature of innate immune responses in the pathogenesis of retinal neovascularization.² Indeed, all three crucial immune response components, including immune cells, mediators, and the complement system, contribute to disease progression. Mononuclear phagocytes, among other immune cells including circulating blood monocytes, tissue-resident macrophages, dendritic cells, and microglia, have been indicated to play a central role in immune responses of healthy and diseased retinae.⁶ In the event of an aggravated insult to the retina, such as in PDR, microglia activation persists and accumulates at the site of damage, failing to restore the homeostatic state of the retina. This results in consistent inflammation and thus induces the secretion of chemokines such as C-C motif chemokine ligand 2 (CCL2), which attracts infiltration of monocytes and macrophages from the leaked blood–retinal barrier and the newly formed blood vessels.^7,8 Moreover, elevated levels of inflammatory cytokines, such as TNF-α, IL-8, CCL2, IL-1β, and IL-6, have been found in the vitreous fluid of patients with PDR,^9–11 suggesting increased activity of phagocytes and macrophages. 

In this study, we investigated the transcriptomic features and cellular landscape of the immune response in the retina of a rat model with ischemic retinopathy. We found that the ischemic retina is transcriptionally characterized by increased immunoregulatory interactions, interleukin-relevant signaling, and complement cascade. Through deconvolution analysis, we demonstrated that myeloid cells, specifically M2 macrophages, majorly infiltrate the ischemic retina. Functional enrichment analysis of differentially expressed genes (DEGs) further revealed that the dysregulated genes are associated with leukocyte-mediated immunity, regulation of cytokine production, myeloid leukocyte activation, and mononuclear cell proliferation, indicating the potential role of M2 macrophages in immune infiltration during retinal neovascularization. Further downstream protein–protein interaction analysis revealed a number of potential hub genes that contribute to M2 macrophage infiltration. Deconvolution analysis of transcriptomic data also revealed increased M2 macrophage levels in the macula of patients with both PDR and diabetic macular edema (DME). This study illuminates a useful method for using transcriptomic data to deconvolute not only immune cells but other cell types from bulk ocular tissues, potentially providing a new way to study the complicated pathogenesis of many ocular diseases. 

Materials and Methods

Animals

Pregnant Sprague Dawley rats (female, 12–14 weeks old) were supplied by the Animal Resources Centre, Perth, Australia, and were housed in standard cages with free access to food and water under a cycle of 12 hours light (50 lux illumination) and 12 hours dark (<10 lux illumination) in a temperature-controlled environment to minimize possible light-induced damage to the eye. The protocol used in the study was approved by the Animal Ethics Committees of St. Vincent's Hospital (004/16). All animal studies were performed according to the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research. 

Rat Model of Oxygen-Induced Retinopathy

The oxygen-induced retinopathy (OIR) model was induced in rats as previously described.¹² Briefly, newborn Sprague Dawley rats (within 12 hours of birth; postnatal day 0 [P0]) and their nursing mothers were subjected to daily cycles of 80% O₂ for 21 hours and room air for 3 hours in a custom-built and humidity-controlled (55%–75%) chamber until P14. Animals were returned to room air until P18, when they were sacrificed to harvest retinae. An oxygen controller (ProOx 110; BioSpherix, Parish, NY, USA) was used to monitor and regulate oxygen levels. Rat retinae were dissected and stained with Life Technologies Isolectin GS-IB4–Alexa Fluor 488 (5 µg/mL, I21411; Thermo Fisher Scientific Australia, Scoresby, VIC, Australia) to confirm the induction of OIR before processing for bulk RNA sequencing (RNA-seq). 

Extraction of Total RNA

Total retinal RNA from OIR rats (n = 3) and control rats (n = 3) was isolated using the RNeasy Mini Kit (74104; QIAGEN, Hilden, Germany). The quantity of RNA was measured using spectrophotometry, and the integrity of total RNA was confirmed using a bioanalyzer at the Australian Genome Research Facility (Melbourne, VIC, Australia). The quality of RNA was determined by RNA integrity number (RIN). A RIN of >7 and concentration of >20 ng/µL were considered acceptable for RNA sequencing. 

Bulk RNA-Seq Analysis

Analysis was undertaken in triplicate of total RNAs (1 µg) from three control and three OIR rats at P18, prepared as per the manufacturer’s instructions. cDNA libraries were sequenced by using the HiSeq 2500 RNA-seq platform (Illumina, San Diego, CA, USA) as 100-base pair single-end chemistry at the Australian Genome Research Facility. The raw fastq files of single-end reads were removed from adapter sequences and the low-quality reads were dropped by TrimGalore 0.6.7 (Babraham Bioinformatics, Cambridge, UK). Filtering parameters were set as follows: -q 25 –length 50 -e 0.1 –stringency 5. The trimmed reads were subjected to alignment with default settings using HISAT2-2.2.1. The reference genome used for rats was rn6. Aligned RNA-seq data were counted over gene exons using featureCounts 2.0.1. Genes were annotated as per the Rattus_norvegicus.Rnor_6.0.104 annotation file. 

Data Accession

The RNA-seq data for the retinae of control rats and OIR rats have been deposited in the Gene Expression Omnibus (GEO) under GSE104588. To access gene expression in the macula of patients with both PDR and DME, we obtained the RNA-seq read count from GSE160306 (accessed October 25, 2021). We specifically compared macular samples from these patients (GSM4871170, GSM4871171, and GSM4871175) with macular samples from healthy controls (GSM4871112, GSM4871123, GSM4871134, GSM4871145, GSM4871155, GSM4871166, GSM4871177, GSM4871188, GSM4871190, and GSM4871113). Supporting data are available in the Supplementary Material files and from the corresponding authors upon request. 

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) was run on both the rodent and human transcriptome dataset against canonical pathway gene set collections (c2.cp.v7.2; Broad Institute of MIT and Harvard, Cambridge, MA, USA) by using 1000 gene set permutations. Only gene sets with a false discovery rate < 0.05 were considered to be significantly enriched. The raw counts were loaded into R 4.1.0 (R Foundation for Statistical Computing, Vienna, Austria) for statistical analysis and the pheatmap function was used to perform hierarchical clustering analyses and for data representation. 

Evaluation of Immune Cell Infiltration by CIBERSORTx

Normalized gene expression data from bulk RNA-seq were used to infer the estimated proportions of infiltrating immune cells using the CIBERSORTx algorithm, a deconvolution algorithm that uses a set of reference gene expression values.¹³ We used the signature matrix LM22, defining 22 infiltrating immune cells, including subsets of T cells, B cells, monocytes, macrophages, dendritic cells, eosinophils, and neutrophils, among others. The raw counts were transformed to transcripts per million, similar to those resulting from microarrays. Normalized data were uploaded to the CIBERSORTx web portal (https://cibersortx.stanford.edu/), and the algorithm was run using the LM22 signature matrix at 1000 permutations in B-mode, using batch correction to adjust the gene expression of the bulk RNA-seq.¹⁴ 

Identification of DEGs and Functional Correlation Analysis

For visualization of adjusted expression levels, the raw counts were normalized and plotted with the DESeq2 package in R 4.1.0. The effect of intersample correction was demonstrated using a two-dimensional principal component analysis (PCA) cluster plot. DEG analyses were performed by DESeq2,¹⁵ and a volcano plot of DEGs was drawn using the EnhancedVolcano package. DEGs were determined at a significance threshold of adjusted P < 0.05 and |log₂FC| > 1. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the clusterProfiler package.¹⁶ 

Protein–Protein Interaction Network Analysis

The DEGs were uploaded to the Search Tool for the Retrieval of Interacting Genes (STRING) database (version 11.0b) to obtain the protein–protein interaction network.¹⁷ The resulting outcome was imported into Cytoscape 3.8.2, followed by cluster analyses of DEGs using the Molecular Complex Detection (MCODE) plug-in, an algorithm to find highly interconnected regions in a network.¹⁸ Genes in the cluster with MCODE score greater than 10 were selected for downstream analysis, including GO and KEGG enrichment analysis. The Cytoscape cytoHubba plug-in was further used to calculate the degree of interaction of selected genes.¹⁹ Hub genes was determined from the top 10 genes identified by cytoHubba. 

Statistical Information Summary

Comprehensive information on the statistical analyses used is included throughout this text, including the figure legends and results. The methods, significance, and P values are described. If not specified, P < 0.05 was considered statistically significant. 

Results

Active Immune Responses in the Retinae of OIR Rats

To profile the transcriptome landscape of retinal neovascularization by bulk RNA-seq, we used a rat model of OIR that yields ischemic avascular zones and preretinal neovascularization resembling clinical observations (Fig. 1A). Significant oxygen-level changes cause a relatively hypoxic condition, resulting in maximal preretinal neovascularization occurring between vascular and non-vascular zones in OIR rats at P18 (Fig. 1B). To gain further insights into the key molecular signaling or biological pathways involved in the pathological process of retinal neovascularization, the transcriptome data were interrogated with GSEA, a knowledge-based approach for interpreting genome-wide expression profiles.²⁰ Surprisingly, all of the top 10 gene sets (Reactome pathway database) positively enriched by GSEA were associated with immune responses, such as immune cell interactions, interleukin or interferon signaling, and complement cascade, suggesting the presence of rather strong immune responses in the retinae of OIR rats (Fig. 1C). Moreover, all of the top 10 gene sets negatively enriched were closely related to the neuronal system and signal transmission, indicating severely aberrant retinal function in OIR rats. Overall, our data imply that immune responses might play a critical role in retinal neovascularization and suggest that various immune cell types are involved in this pathological process. 

Figure 1.

 

Activation of immune response in retinal neovascularization. (A) Schematic diagram of the Sprague Dawley rat model of OIR and workflow of RNA-seq data analysis. Created with Biorender.com. (B) Representative images of retinal whole mounts of OIR and control rats. Vasculature was visualized by IB4 labeling (green). Scale bars: 1 mm (left panel) and 100 µm (right panel). (C) GSEA revealed that immunoregulatory interactions, interleukin-relevant signaling, and complement cascade were positively enriched, whereas neuronal systems and signal transmission were negatively enriched in OIR rats.

M2 Macrophages Are the Major Immune Cells That Infiltrate the Retinae of OIR Rats

In order to understand different immune cell types and to assess their expression in the retinae of OIR rats, we performed immune cell infiltration analysis by CIBERSORTx, a machine learning method that enables inference of cell-type–specific gene expression profiles from bulk tissue transcriptome data without physical cell isolation.¹³ The estimated proportions of immune cells in both OIR and control retinae are shown in the bar plot and heatmap (Fig. 2A), which indicate differential immune infiltrates between the two groups. We used the Kruskal–Wallis test to compare the fraction of each immune cell type between OIR rats and control rats and found that M2 macrophages, eosinophils, plasma cells, and CD8 T cells significantly infiltrated the retinae in OIR rats, among 22 immune cells that were deconvoluted by CIBERSORTx (Fig. 2B). Notably, we observed an evidently increased proportion of resting CD4 memory T cells in controls compared with OIR rats. Further analysis with Welch's t-test suggested that only M2 macrophages were significantly increased in OIR rats compared with control rats (Fig. 2C), whereas no significant differences for other cell types were observed between the two groups (Supplementary Fig. S1). To identify co-expression patterns among significantly altered immune cells, we performed Spearman correlation analyses to evaluate possible relationships. The results showed that M2 macrophages had a strong positive correlation with eosinophils (correlation coefficient r = 0.72) (Fig. 2D), whereas both M2 macrophages and eosinophils had a negative correlation with resting CD4 memory T cells (r = −0.77 for M2 macrophages and r = −0.88 for eosinophils), suggesting that the innate immune response (M2 macrophages and eosinophils) rather than adaptive immune response (resting CD4 memory T cells) is more prominent in OIR rats. 

Figure 2.

 

The immune cellular landscape of the retinae of OIR rats. (A) Estimated proportions of each of the 22 immune cells in the retinae of OIR and control rats as determined by CIBERSORTx. (B) Comparison of estimated proportion of each immune cell type between OIR rats (n = 3) and control rats (n = 3). Kruskal–Wallis tests showed that M2 macrophages were significantly increased in OIR rats. *P < 0.05, ns, no significance; ?, cells were not detected so comparison was not possible. (C) Welch's t-tests showed that the M2 macrophages were significantly increased in OIR rats. (D) Correlation heatmap of 22 immune cell types. Values greater than 0.3 indicate that the two cell types are positively correlated, whereas values less than −0.3 indicate that the two cell types are negatively correlated. NA, not available.

Functional Enrichment Analysis of DEGs Between OIR and Control Rats

Having observed the differential abundance of M2 macrophages in OIR versus control rats using CIBERSORTx, we used transcriptomic data to further validate these estimates. PCA can be used to assess the consistency of biological repetition and the difference between groups. Our PCA results demonstrated a clear difference between the OIR rats and control rats (Fig. 3A). We identified 541 DEGs between the OIR rats and control rats (adjusted P < 0.05; |log₂FC| >1) (Supplementary Table S1), including 444 upregulated genes and 97 downregulated genes (Fig. 3B). GO enrichment analysis revealed that DEGs were primarily associated with immune responses in the biological process, including leukocyte-mediated immunity, positive regulation of cytokine production, leukocyte migration and proliferation, and activation of immune responses (Fig. 3C). Of note, GO terms in association with myeloid cells, such as myeloid leukocyte activation and mononuclear cell proliferation, were ranked among the top 10. Interestingly, KEGG enrichment analysis also demonstrated that a number of top pathways were enriched in the function of myeloid cells, such as osteoclast differentiation, phagosome, neutrophil extracellular trap formation, and chemokine signaling pathway (Fig. 3D). In addition to our immune infiltration data, these DEG enrichment results further suggest that myeloid cells, particularly M2 macrophages, play a principal role in immune responses during retinal neovascularization. 

Figure 3.

 

Functional enrichment analysis of DEGs. (A) PCA suggested that the transcriptomic profiles of OIR rats and controls are evidently distant. (B) Volcano plot of 541 DEGs between OIR rats and control rats, including 444 upregulated genes and 97 downregulated genes. (C) GO and (D) KEGG enrichment analysis of DEGs. BP, biological process; CC, cellular component; MF, molecular function.

Identification of Hub Genes Involved in M2 Macrophages Infiltration in OIR

To identify the hub genes involved in M2 macrophages infiltration in OIR, a protein–protein interaction network of DEGs was constructed using the STRING online database. We identified 452 nodes and 3933 edges among the DEGs. The STRING data were imported to Cytoscape and processed with MCODE, resulting in 12 protein clusters (Fig. 4A). We selected the top three clusters with MCODE score above 10 (visualized with Cytoscape) (Supplementary Fig. S2) for further GO enrichment analysis. The results showed that genes involved in the selected clusters were mainly associated with leukocytes, particularly myeloid leukocyte activities, in biological processes such as cell migration, activation and chemotaxis (Fig. 4B). We subsequently selected all GO terms relevant to macrophages, resulting in 26 unique and associated genes (Table). To further pinpoint hub genes most relevant to macrophages, we applied the Cytoscape plugin cytoHubba to rank the genes according to the degree of their proteins, a topological analysis method to assess the essentiality of genes.¹⁹ This analysis resulted in 10 top-scored genes (Fig. 4C), including Ccl2, integrin alpha M (Itgam), toll-like receptor 2 (Tlr2), vascular endothelial growth factor A (Vegfa), colony-stimulating factor 1 receptor (Csf1r), C-X-C motif chemokine ligand 12 (Cxcl12), colony stimulating factor 1 (Csf1), allograft inflammatory factor 1 (Aif1), triggering receptor expressed on myeloid cells 2 (Trem2), and C-C motif chemokine ligand 3 (Ccl3). A heatmap further demonstrated that expression of these genes (except for Cxcl12) was evidently higher in OIR rats compared with control rats (Fig. 4D), suggesting that they are potentially crucial for the activity of M2 macrophages in retinal neovascularization. 

Figure 4.

 

Identification of hub genes associated with infiltration of M2 macrophages in retinal neovascularization. (A) Gene cluster identified by MCODE with DEGs. The top three clusters (score > 10) were selected for downstream analysis. (B) GO analysis of all genes from the selected three clusters. Genes from GO terms related to the macrophages identified were selected to identify hub genes by the Cytoscape cytoHubba plug-in. (C) Hub genes identified by cytoHubba potentially contribute to infiltration of M2 macrophages. (D) Heatmap of expression of the identified hub genes.

Table.

 

GO Terms Relevant to Macrophages and Their Corresponding Genes, Analyzed With the Top Three Protein Clusters by STRING

 

M2 Macrophages Are the Major Immune Cells to Infiltrate in the Macula of Patients With Both PDR and DME

To understand the immune responses and immune cells involved in patients with retinal neovascularization, we analyzed a published RNA-seq dataset from the macula of patients with both PDR and DME. GSEA analysis revealed a number of positively enriched gene sets in the macula of patients with both PDR and DME, including complement cascade, immune cell interactions, and interleukin signaling, whereas gene sets negatively enriched were closely related to the neuronal system and signal transmission (Fig. 5A). We further estimated the proportion of immune cells in these human samples by CIBERSORTx as shown in Figure 5B. Among 22 immune cells, M2 macrophages and memory B cells were evidently increased and eosinophils were decreased in patients with both PDR and DME compared with healthy controls, as analyzed by the Kruskal–Wallis test (Fig. 5C). Further analysis with Welch's t-test indicated that only M2 macrophages were significantly increased in the macula of the patients (Fig. 5D). Moreover, we also observed that expression of the 10 hub genes identified in OIR rats was clearly increased in the macula of these patients (Fig. 5E). Together, these results generated from human samples showed a close similarity to those from OIR rats, confirming that active M2 macrophages are the major immune cells involved in the pathogenesis of retinal neovascularization. 

Figure 5.

 

The immune cellular landscape of the macula of patients with PDR and DME. (A) GSEA revealed that the initial triggering of complement, complement cascade, immunoregulatory interactions, and interleukin-relevant signaling was positively enriched, whereas the neuronal system and signal transmission were negatively enriched in the macula of patients with PDR and DME. (B) Estimated proportions of each of the 22 immune cells in the macula of patients with both PDR and DME (n = 10) compared with healthy controls (n = 3) by CIBERSORTx. (C) Comparison of the estimated proportion of each immune cell type between patients with PDR and DME and healthy controls. Kruskal–Wallis tests showed that M2 macrophages and memory B cells were significantly increased in the macula of patients with PDR and DME, whereas eosinophils were decreased. *P < 0.05; ns, no significance; ?, cells were not detected so comparison was not possible. (D) Welch's t-tests revealed significantly increases in M2 macrophages in the macula of patients with PDR and DME. (E) Heatmap of expression of the hub genes (identified in OIR rats) in the macula of patients with PDR and DME. Expression of most of the genes in the human samples shared profiles with those in OIR rats.

Discussion

Here, we systemically analyzed the transcriptome and immune infiltration of the retinae of rats subjected to ischemia-induced retinopathy. GSEA results demonstrated the strong immune responses in relation to retinal neovascularization, including immunoregulatory interactions, interleukin-relevant signaling, and the complement cascade. Using CIBERSORTx to estimate the proportions of 22 immune cells in bulk tissues, we revealed that M2 macrophages were significantly increased in the OIR rat models. In line with GSEA results, functional enrichment analysis of DEGs further pinpointed that leukocyte-mediated immunity and cytokine production was involved in retinal neovascularization. Protein–protein interaction analysis resulted in three major protein clusters, and functional enrichment analysis further suggested that myeloid leukocyte-mediated immune responses play a key role. Macrophage-relevant genes from the selected protein clusters were further screened to identify a number of hub genes that contribute to M2 macrophage infiltration, including Ccl2, Itgam, Tlr2, Vegfa, Csf1r, Cxcl12, Csf1, Aif1, Trem2, and Ccl3. Moreover, analysis of transcriptomic data from the macula of patients with both PDR and DME also suggested that M2 macrophages are the major immune cells involved in the pathogenesis of retinal neovascularization. 

It has been suggested that the innate immune system is an early responder to environmental perturbations to maintain retinal homeostasis and visual function.²¹ Retinal homeostasis is retained majorly through the blood–retinal barrier and several specialized cell types, such as microglia and macrophages. In the event of an insult, metabolic dysregulation activates microglia which in turn enhances the innate immune response by releasing inflammatory mediators, including IL-1β and TNF-α and recruits other immune cells to eliminate the inciting stressor.^22–24 However, this microglia response becomes less effective with persistent perturbation, such as diabetic retinopathy, and eventually contributes to chronic inflammation. Macrophages are one of the innate immune cells suggested to participate in the regulation of the unique retinal immune responses, playing a role as a pseudo-barrier in clearing proteins and debris.² However, it remains unclear which subsets of macrophages are involved and how they contribute to the retinal immune responses. Through analysis of the retinal transcriptome from a rat model of ischemic retinopathy, we demonstrated that the retina undergoes severe immune responses, mainly mediated by interleukin signaling, interferon signaling, and complement cascade. Consistent with these findings, our analysis of the published human retinal transcriptome (macula from patients with both PDR and DME) revealed that gene expression in the macula is closely related to the immune response-related signaling found in the rodent model, indicating that strong immunogenicity occurs during retinal neovascularization. 

It is known that macrophages play an important role not only in immune response but also at each step of the angiogenic cascade, from forming neovascular sprouting to the remodeling of the vascular plexus and vessel maturation.²⁵ Macrophages are heterogeneous and characterized differently by their function and activation route. Under pathophysiological conditions, such as inflammation and angiogenesis, macrophages are activated and can be classically categorized as M1 and M2 phenotypes.²⁶ M1 macrophages that are induced by IFN-γ, TNF-α, and bacterial lipopolysaccharides, are characterized as proinflammatory, whereas M2 macrophages that can be stimulated by IL-4, IL-10, and IL-13 are characterized as antiinflammatory.^27,28 Moreover, the complement system consisting of proinflammatory mediators, such as C3a and C5a, has been indicated to induce macrophage infiltration in angiogenesis and chronic inflammation.²⁹ Indeed, studies have shown that M2 but not M1 macrophages express higher angiogenic growth factors and cytokines that promote angiogenesis both in vitro and in vivo.^30,31 Nevertheless, the role of macrophages in retinal neovascularization remains unclear. Interestingly, our GSEA results showed positive enrichment of IL-4, IL-10, and IL-13 signaling and complement cascade in both OIR rats and patients with PDR, indirectly suggesting a potential role of M2 macrophages in retinal neovascularization. This assumption was confirmed by the analysis with CIBERSORTx, which identified a significantly increased proportion of M2 macrophages in the retinae of OIR rats. This finding is highly consistent with results from a previous report that demonstrated that predominant M2 macrophages promote retinal pathological neovascularization by producing secreted factors in a mouse model of OIR.³² A study also reported increased M2 macrophages in the vitreous humor of patients with PDR compared with non-diabetic controls, induced by macrophage colony-stimulating factor (M-CSF) and IL-13.³³ In the present study, we used CIBERSORTx to demonstrate for the first time, to the best of our knowledge, that M2 macrophages are also significantly increased in the macula of PDR patients. Although we observed increased numbers of a few other immune cell types in OIR rats, such as eosinophils, plasma cells, and CD8 T cells, there was no significant difference between the OIR rats and the control rats. Such changes of these immune cell types were also not observed in human PDR samples, which is likely due to the limited number of samples and increased complexity of PDR and DME pathogenesis in humans. Nevertheless, these results clearly indicate that M2 macrophages are crucial to the innate immune response in the pathogenesis of retinal neovascularization. 

Our results show that Ccl2, Itgam, and Tlr2 are the top three hub genes with the highest degree of interaction, suggesting their potential role in M2 macrophage infiltration in retinal neovascularization. Not only did we observe increased expression levels of these genes in the rodent model of retinal neovascularization, but we also identified a similar phenomenon in the human transcriptome. Ccl2 has been reported as an important chemokine secreted by activated microglia and monocytes from the vascular supply, which recruits monocytes and macrophages to take part in the pathogenesis of retinal degeneration and choroidal neovascularization. In the age- and light-induced atrophic age-related macular degeneration mice models (Cx3cr1^−/−) that feature subretinal inflammation and photoreceptor degeneration, ablation of Ccl2 or its receptor reduces monocyte infiltration and retinal degeneration, suggesting its regulatory role in the innate immune response of insulted retina.³⁴ Moreover, a study found that CCL2 protein levels were increased significantly in OIR mice, resulting in inflammation of the ischemic retina, possibly by modulation or recruitment of macrophages.³⁵ Itgam is essential to macrophage and monocyte activities, such as adhesive interactions and macrophage polarization.^36,37 It may enable antiinflammatory macrophage polarization via activation of transcription factors such as Stat3 to promote angiogenesis.³⁶ In fact, increased protein levels of Itgam were observed in retinal neovascular tufts in Vldlr^–/– mice, characteristic of abnormal subretinal neovascularization. It has been suggested that Tlr2 secreted by macrophages regulates the expression of angiogenic factors, such as VEGF, and may be important in ischemia when TLR agonists are present.³⁸ A recent study reported that Tlr2 knockout attenuated retinal neovascularization in OIR mice by decreasing mRNA expression of VEGF and a number of inflammatory factors, including TGF-β1, b-fibroblast growth factor (bFGF), and IL-6, suggesting that inflammatory responses induced by Tlr2 are associated with angiogenesis. Other hub genes including Vegfa, Csf1r, Cxcl12, Csf1, Aif1, Trem2, and Ccl3 have been shown to play a role in the regulation of recruitment, growth, polarization, and differentiation of M2 macrophages.^39–44 However, the role of each hub gene predicted in this study in the regulation of M2 macrophages and how they contribute to the pathogenesis of retinal neovascularization remain unclear. 

There are some limitations to this study. First, the number of samples in each group for RNA-seq was relatively small. Greater sample numbers are warranted to verify our results. Second, this study was based purely on bioinformatic analyses using various algorithmic tools, for which biological validation will be necessary, particularly for the predicted hub genes. In addition, further studies should be performed to assess if the cellular phenotypic change results in the altered expression of these hub genes. 

Conclusions

We found that M2 macrophages are the principal immune cells involved in the innate immune response of retinal neovascularization using the gene expression-based deconvolution method. We also predicted several hub genes associated with the identified immune cell types through various algorithmic tools. Our results offer unforeseen insights into the landscape of immune cells associated with retinal neovascularization and provides a new and accessible approach to investigate the role of immune cells in many other ocular disorders, through analysis of transcriptomic data. 

Acknowledgments

Supported by grants from the National Health and Medical Research Council of Australia (1185600). The Centre for Eye Research Australia receives Operational Infrastructure Support from the Victorian Government. All datasets generated for this study are freely available in the Supplementary Materials or from the corresponding author upon request. 

Disclosure: J.-H. Wang, None; S. Kumar, None; G.-S. Liu, None