Retinal Transcriptome and Cellular Landscape in Relation to the Progression of Diabetic Retinopathy

Abstract

Purpose: Previous studies that identify putative genes associated with diabetic retinopathy are only focusing on specific clinical stages, thus resulting genes are not necessarily reflective of disease progression. This study identified genes associated with the severity level of diabetic retinopathy using the likelihood-ratio test (LRT) and ordinal logistic regression (OLR) model, as well as to profile immune and retinal cell landscape in progressive diabetic retinopathy using a machine learning deconvolution approach.

Methods: This study used a published transcriptomic dataset (GSE160306) from macular regions of donors with different degrees of diabetic retinopathy (10 healthy controls, 10 cases of diabetes, 9 cases of nonproliferative diabetic retinopathy, and 10 cases of proliferative diabetic retinopathy or combined with diabetic macular edema). LRT and OLR models were applied to identify severity-associated genes. In addition, CIBERSORTx was used to estimate proportional changes of immune and retinal cells in progressive diabetic retinopathy.

Results: By controlling for gender and age using LRT and OLR, 50 genes were identified to be significantly increased in expression with the severity of diabetic retinopathy. Functional enrichment analyses suggested these severity-associated genes are related to inflammation and immune responses. CCND1 and FCGR2B are further identified as key regulators to interact with many other severity-associated genes and are crucial to inflammation. Deconvolution analyses demonstrated that the proportions of memory B cells, M2 macrophages, and Müller glia were significantly increased with the progression of diabetic retinopathy.

Conclusions: These findings demonstrate that deep analyses of transcriptomic data can advance our understanding of progressive ocular diseases, such as diabetic retinopathy, by applying LRT and OLR models as well as bulk gene expression deconvolution.

Diabetic retinopathy (DR) is a progressive retinal complication of diabetes, causing significant visual impairment on a global scale.¹ Microvascular lesions of the retina have been used to evaluate and classify the severity level of DR. There are two broad categories, including the early stage of nonproliferative diabetic retinopathy (NPDR) and the advanced stage of proliferative diabetic retinopathy (PDR).¹ NPDR is manifested by increased vascular permeability and capillary occlusion, whereas PDR is characterized by retinal neovascularization in addition to NPDR features.² An important additional category of DR is diabetic macular edema (DME), the most common cause of vision loss in patients with DR.² Genetic studies, including candidate gene studies, linkage studies, powerful genomewide association studies (GWAS), and meta-analysis, have identified putative candidate gene or single nucleotide polymorphisms associated with DR.^3,4 Although these studies provide critical information that helps to interpret the mechanisms of DR, they mostly study one specific type of DR, thus the resulting putative genes are not likely to be associated with DR progression. 

In this study, we assess the macular transcriptomic data from donors with different severity levels of DR and define the genes changed with DR progression. The landscape of immune cells and retinal cells in relation to DR progression has also been established by using a machine learning deconvolution algorithm, CIBERSORTx. 

Materials and Methods

Data Accession

RNA-Seq raw counts were obtained from Gene Expression Omnibus (GSE160306, accessed on March 28, 2022). Raw counts were converted to transcripts per million (TPM) in R (version 4.1.0) for downstream analyses. Samples from the human macular region were selected from the raw dataset using the subset function from the R package Seurat (version 4.0.5).⁵ Detailed clinical information, including patient gender, age, and disease stages, were accessed from GSE160310. Single-cell transcriptomes of the human neural retina were directly downloaded from Human Cell Atlas (HCA; https://www.humancellatlas.org/, E-MTAB-7316, accessed on January 15, 2022).⁶ 

Identification of Severity-Associated Genes in the Macular Region From Donors With DR at a Different Stage

To identify genes related to the progressive severity of DR, the selected RNA-Seq raw count was analyzed by DESeq2 (version 1.32.0),⁷ using the likelihood-ratio test (LRT), as we previously described.⁸ Briefly, gene expression was determined by three variables due to limited access to clinical labeling, consisting of gender, age, and disease status (severity). To simplify the design in the LRT, the progressive disease status was indexed using integer scaling: control-1, diabetic-2, NPDR-3, and NPDR/PDR+DME-4. Gender and age as covariants were regressed out in the LRT model, leaving disease status as the main factor to affect gene expression. The batch correction was performed using the removeBatchEffect function in limma (version 3.48.3).⁹ Significant severity-associated genes were determined by adjusted P < 0.05, controlled for gender and age. The resulting genes were then scaled to z-score and clustered using the degPatterns function from the R package DEGreport (version 1.28.0). Gene clusters with consistent and progressive changes were selected for further analyses. 

To further screen significant severity-associated genes, the TPM values of the genes from the selected cluster were analyzed by R package MASS (version 7.3-54) and ordinal (version 2019.12-10) using the ordinal logistic regression (OLR) model.¹⁰ Using Cumulative Link Model function in OLR,¹⁰ we controlled for gender and age and identified input genes that were associated with disease severity. These severity-associated genes were further determined at P < 0.05. The resulting genes were represented in a forest plot and ranked according to their beta coefficient. A heatmap of severity-associated genes was generated by the heatmap function from the R package NMF (version 0.24.0).¹¹ 

To assess the gene-gene correlation of severity-associated genes, the raw counts were processed by the Variance Stabilizing Transformation function in DESeq2.⁷ followed by statistical adjustment for gender and age using the Remove Batch Effect function in limma (version 3.48.3).⁹ The correlation heatmap was plotted using R package corrplot (version 0.92).¹² 

Functional Enrichment Analysis and Gene Expression Visualization

Functional enrichment analyses (Gene Ontology Biological Processes) for significant severity-associated genes were performed with Metascape (https://metascape.org/gp/index.html#/main/step1),¹³ which provides more frequently updated bioinformatics analyses than DAVID.¹⁴ 

To visualize the expression of selected genes, the raw counts were processed by the Variance Stabilizing Transformation function in DESeq2.⁷ followed by statistical adjustment for gender and age using the Remove Batch Effect function in limma (version 3.48.3).⁹ Tukey boxplots were used to compare gene expression from patients with different DR statuses, with interquartile range boxes and 1.5 times interquartile range whiskers. The 2-tailed Mann-Whitney test was applied to assess pairwise comparisons in the plots, whereas the Kruskal-Wallis test was used to perform statistical comparison across all groups of different disease statuses. 

Deconvolution Analyses of the Immune Cells and Retinal Cells From Donors With DR

The CIBERSORTx, a machine learning deconvolution algorithm that enables inference of cell-type-specific gene expression profiles,¹⁵ was applied to infer the estimated proportions of infiltrating immune cells or retinal cells of each bulk retinal transcriptome. We used the LM22 signature matrix to define 22 infiltrating immune cells.¹⁶ Normalized data (TPM value) of the transcriptome were uploaded to the CIBERSORTx web portal (https://cibersortx.stanford.edu/), with the algorithm run using LM22 signature matrix at 100 permutations with B-mode batch correction. 

To infer the retinal cellular profile, we used a previously derived retinal signature matrix using the HCA scRNA-Seq dataset of the human retina.⁸ In brief, the retinal signature matrix contains rod, Müller glia, bipolar, cone, amacrine, microglia, ganglion, astrocytes, and horizontal cells. The retinal signature matrix was used to impute retinal cellular fraction from the bulk transcriptome at 100 permutations with S-mode batch correction. Statistical significance of cellular numbers associated with severity was assessed by an ordinal logistic regression model using the nonparametric Kruskal-Wallis test that allows for multifactorial designs. The estimated coefficients and P value were presented as a forest plot with 95% confidence intervals. The estimated proportions of cell types were adjusted for gender and age to assess the effects by disease severity by extracting the residuals from fitting a generalized linear model with variables. 

Statistical Analyses

Detailed information of statistical analyses was described in various places, including results and figure legends, where the test methods, significance, P values, error bars, and coefficients were included. 

Results

Identification of Severity-Associated Genes in the Macula of Patients With DR at Different Clinical Stages

We analyzed a publicly available transcriptional dataset of the macular region from the retina of patients at different stages of DR and healthy controls. A total of 39 macular RNA-Seq profiles from healthy controls (n = 10) and cases from early to late DR stage (n = 29) were analyzed (Fig. 1A). The patient groups were previously classified according to the severity of DR progression, including diabetic (n = 10), NPDR (n = 9), NPDR with DME (n = 7), and PDR with DME (n = 3).¹⁷ Of which, patients with NPDR with DME or PDR with DME were combined as one group of NPDR/PDR+DME (n = 10), as the previous study described.¹⁷ 

Figure 1.

 

Identification of severity-associated genes in the macula from patients with degrees of DR. (A) Demographics of the human macular RNA-Seq profiles, detailed by age and different degrees of DR. (B) Tukey boxplots (interquartile range [IQR] boxes with 1.5 × IQR whiskers) of severity-associated genes found in the macula of patients with DR. The expression level of severity-up genes increases with the severity of DR (left, n = 63). Severity-up genes were defined with DESeq2 two-sided likelihood-ratio test (adjusted P < 0.05) by controlling for gender and age. Adjusted gene expressions were shown as z-score. (C) Forest plot of severity-associated genes significantly correlated with DR severity (P < 0.05, n = 50). Positive coefficients indicate severity-associated genes increase in expression with DR severity. Statistical analyses of the correlation between the severity-associated genes and the DR progression were assessed by a nonparametric ordinal logistic regression model that controls for gender and age. Point sizes are scaled by statistical significance. Error bars represent 95% confidence intervals.

Gene expression levels in the retina were largely affected by age, gender, and disease stages.¹⁸ By controlling for age and gender, we identified a cluster of genes (n = 63) related to DR progression from the early to late stage (Fig. 1B) using an LRT.⁷ Of note, the expression level of these genes was increased with the severity of DR. We therefore defined them as severity-up genes. To further screen genes significantly changed with the severity of DR, we assessed the correlation between severity-up genes and the severity of DR by using the OLR model that controls for gender and age. Our results identified 50 genes progressively increasing their expression level with the severity of DR (Fig. 1C, hereafter referred to as severity-associated genes, P < 0.05). The remaining 13 severity-up genes were not significantly correlated to the severity of DR, thus excluding from the downstream analyses. A heatmap clearly demonstrated that the expression levels of 50 severity-associated genes were progressively increased with the severity of DR from the early (diabetic) to late stage (NPDR/PDR+DME), after adjustment for gender and sex (Fig. 2A). Correlation analysis among these genes in terms of expression level revealed a cluster of genes (n = 25) having the most interactions and strong correlations (|r| > 0.6) with other genes, suggesting their potentially important roles in the progression of DR (Fig. 2B). 

Figure 2.

 

Expression profile and correlation analyses of severity-associated genes. (A) A heatmap of the expression level of 50 severity-associated genes in all samples across different groups, adjusted for gender and age. (B) Heatmap of the correlation matrix across 50 severity-associated genes. Pearson's correlation was calculated among the severity-associated genes to show the co-expression patterns of genes in the heatmap. Color key denotes the Pearson's correlation coefficient between genes. There is one clear cluster (highlighted in red and black) of severity-associated genes with a strong positive correlation, adjusted by gender and age.

Functional Enrichment Analyses of the Severity-Associated Genes and Identification of Key Genes Involved in the Progression of DR

To determine the functional roles of the severity-associated genes in the pathogenesis of DR, functional enrichment analysis (Gene Ontology [GO] biological processes) was performed with Metascape. The results demonstrated significant enrichments mostly related to inflammation and immune response. Of which, negative regulation of cell differentiation, response to amyloid-beta, and inflammatory response rank top three among other GO terms (Fig. 3A). Genes involved in the top three GO terms were combined to identify key genes that contribute to the progression of DR development, resulting in seven genes with significant changes in expression with the severity of DR after controlling for gender and age (Fig. 3B), including cyclin D1 (CCND1), Fc gamma receptor IIb (FCGR2B), matrix metallopeptidase 9 (MMP9), toll like receptor 4 (TLR4), alpha tocopherol transfer protein (TTPA), CKLF like MARVEL transmembrane domain containing 5 (CMTM5), and NLR family pyrin domain containing 1 (NLRP1). To further identify key severity-associated genes involved in the progression of DR, we overlapped the cluster of severity-associated genes (n = 25) identified in Figure 2B with 7 genes identified in the functional enrichment analysis, leading to 3 genes shared by both clusters of genes (Fig. 3C), including CCND1, FCGR2B, and NLRP1. We focused on CCND1 and FCGR2B and investigated their role in immune cells and retinal cells during DR. Given that the expression pattern of NLRP1 was not consistently increased with the severity of DR, we therefore excluded it for further analyses. 

Figure 3.

 

Functional enrichment analyses of severity-associated genes and identification of key genes involved in the DR progression. (A) Functional enrichment analysis (gene ontology annotations of biological processes) of 50 severity-associated genes by Metascape. The top three biological processes highlighted in red were selected for down-stream analysis. Genes with significant changes in these selected annotations include CCND1, FCGR2B, MMP9, TLR4, TTPA, CMTN5, and NLRP1. (B) Tukey boxplots (interquartile range [IQR] boxes with 1.5 × IQR whiskers) showing the expression of CCND1, FCGR2B, MMP9, TLR4, TTPA, CMTN5, and NLRP1. Gene expression values are shown as log-transformed, controlled for gender and age. Statistical significance of difference was assessed by a two-sided Kruskal-Wallis test on the adjusted expression values. (C) Overlapping the clustered genes from Figure 2B (n = 25) and selected genes in Figure 3B (n = 7) resulted in 3 genes, including CCND1, FCGR2B, and NLRP1.

The Immune Cellular Landscape Changed in the Progression of DR

Given that severity-associated genes are strongly correlated to immune response, we next sought to understand the role of immune cells in the retina of patients with DR at different clinical stages. By using CIBERSORTx, we deconvoluted the macular transcriptomic data (Supplementary Fig. S1) with the signature matrix of immune cells, LM22.¹⁹ We further assessed what immune cell types change in proportions with the severity of DR using an OLR model. By controlling for gender and age, we demonstrated that estimated proportions of memory B cells and M2 macrophages were significantly increased with the progression of DR, whereas the estimated proportions of resting mast cells and naïve B cells were significantly decreased with the progression of DR (Fig. 4A, Supplementary Fig. S2). Other immune cell types were neither found in the macular tissue nor their estimated proportions were not significantly altered with the progression of DR. 

Figure 4.

 

The immune cellular landscape of the macula of patients with DR and its relation to severity-associated genes. (A) Forest plot of estimated proportions of immune cells having a significant correlation with the severity of DR (P < 0.05). Positive or negative coefficients indicate the proportion of immune cells increases or decreases with the severity of DR. Statistical analyses of the correlation between the severity of DR and the proportions of immune cells were assessed by a nonparametric ordinal logistic regression model, adjusted for gender and age. Point sizes are scaled by statistical significance. Error bars represent 95% confidence intervals. (B, C) Analyses of the correlation between the expression level of CCND1 or FCGR2B and the estimated proportion of memory B cells or M2 macrophages, respectively, adjusted for gender and age. Pearson's correlation coefficient was used to test the strength of linear relationships between gene expression and the estimated proportion of specific cell types. Grey dots denote individual samples (n = 39). Blue lines denote regression lines.

We further studied the correlation between the expression level of two key genes identified above, CCND1 and FCGR2B, and the proportions of the altered immune cell types, including memory B cells and M2 macrophages, using a Pearson's regression analysis. Our results showed that expression levels of both CCND1 and FCGR2B have a moderate positive (0.4< r <0.6) correlation with the proportions of memory B cells (Fig. 4B) and M2 macrophages (Fig. 4C), respectively. These results suggest that CCND1 and FCGR2B may be involved in the regulation of immune responses during the progression of DR mainly through M2 macrophages. 

The Retinal Cellular Landscape Changed in the Progression of DR

A previous study demonstrated that the retinal structure, particularly the retinal nerve fiber layer, is changed at the early onset of DR.²⁰ However, quantitative changes in retinal cells during the progression of DR have not been reported. Here, we deconvoluted the macular transcriptomes with CIBERSORTx to profile the estimated proportions of each of retinal cell types, using the scRNA-seq data from HCA as a reference (Supplementary Fig. S3). We applied an OLR model to control for gender and age and evaluated if the severity of DR affects the proportions of retinal cell types. The results revealed that the proportions of Müller glia were significantly increased with the progression of DR, whereas the proportions of other retinal cell types were altered but without significant difference (Fig. 5A, Supplementary Fig. S4). To confirm this finding, we further analyzed the expression profile of the severity-associated genes (40 genes were found in the scRNA-Seq dataset) in the human neural retina by using our previous scRNA-Seq data. Of the major retinal cell types, our results showed that most of the severity-associated genes were expressed in a high percentage of Müller glial cells and astrocytes (Fig. 5B). Given the estimated proportion of astrocytes is not significantly changed, we suspect Müller glial cells play an important role in relation to the severity of DR. 

Figure 5.

 

The cellular landscape of the macula of patients with DR and its relation to severity-associated genes. (A) Forest plot of estimated proportions of retinal cells having a significant correlation with the severity of DR (P < 0.05). Positive or negative coefficients indicate the proportion of retinal cells increases or decreases with the severity of DR, respectively. Statistical analyses of the correlation between the severity of DR and the proportions of retinal cells were assessed by a nonparametric ordinal logistic regression model, adjusted for gender and age. Point sizes are scaled by statistical significance. Error bars represent 95% confidence intervals. (B) Heatmap showing the percentage of retina cells expressing severity-associated genes (40 genes were found in the retinal scRNA-Seq dataset from the Human Cell Atlas), scaled by gene across the different cell types. (C) Analyses of the correlation between the expression level of CCND1 or FCGR2B and the estimated proportion of Müller glia, adjusted for gender and age. Pearson's correlation coefficient was used to test the strength of linear relationships between gene expression and the estimated proportion of Müller glia. Grey dots denote individual samples (n = 39). Blue lines denote regression lines.

Similarly, we investigated the relationship between the expression level of two key genes, CCND1 and FCGR2B, and the proportion of Müller glia in DR. Our results showed that the expression level of CCND1 and FCGR2B has a moderate positive (0.4< r <0.6) correlation with the proportion of Müller glia (Fig. 5C), indicating that these genes may be modulated to promote the increase of Müller glia with the severity of DR. 

Discussion

Genetic studies, including linkage analyses, candidate gene association studies, and GWAS, have identified many risk loci or corresponding genes for DR.^3,21,22 However, such studies tended to focus on one specific clinical stage of DR, particularly the advanced stage, PDR, and DME. Therefore, the resulting disease-causing genes from the genetic studies are not necessary to be associated with the progression of DR. Through the application of LRT on transcriptome from patients with different degrees of DR along with OLR modeling, we identified 50 genes that significantly increased in expression level from the early to late stage of DR (severity-associated genes) by controlling for gender and age. Among which, several severity-associated genes have been previously identified to be associated with DR in human or animal models, including TLR4,²³ plasmalemma vesicle-associated protein (PLVAP),²⁴ fibromodulin (Fmod),²⁵ C4A,²⁶ MMP9,²⁷ insulin-like growth factor binding protein 2 (IGFBP2),²⁸ and NLR Family Pyrin Domain Containing 1 (NLRP1),²⁹ whereas none is reportedly increased with the progression of DR. Our results further supported those previous findings, indicating our transcriptome analyses are reliable and comparable to other genetic studies. Many severity-associated genes have not been reported in association with the progression of DR, interestingly, however, they are related to diabetes, such as alpha-2-Macroglobulin Like 1 (A2ML1),³⁰ and C1r,³¹ suggesting their potential role in the pathogenesis of DR. 

Inflammation is a nonspecific response of the immune system to harmful stimuli or stress that includes many functional and molecular mediators.³² Increasing evidence demonstrates that inflammation is crucial to the development of DR. Increased levels of a variety of inflammatory cytokines and chemokines were observed in serum³³ and ocular samples (aqueous and vitreous) from patients with PDR.³⁴ However, it remains unclear how inflammation takes part in the progression of DR due to the difficulty of sample collection. Our functional enrichment analyses revealed that the severity-associated genes are mainly involved in the negative regulation of cell differentiation, inflammation, and immune responses, suggesting these pathological processes are progressive during the development of DR. 

We further identified CCND1 and FCGR2B as two crucial genes involved in the inflammation in DR. CCND1 protein is the critical gatekeeping protein regulating the transition from G1 phase to S phase via the restriction point of the cell cycle.³⁵ A higher mRNA expression level of CCND1 was observed in human diabetic islets by microarray and quantitative real-time PCR (qRT-PCR).³⁶ Heightened hepatic expression of CCND1 was also found in animal models of obesity or diabetes due to hyperinsulinemia,³⁷ suggesting a close relationship between CCND1 and diabetes. The original study with the macular transcriptomic data from patients with DR that we analyzed in the present study also found that CCND1 was one of the genes identified with significant changes in expression levels.¹⁷ Our correlation analysis among severity-associated genes confirmed that CCND1 has a strong interaction with most of other severity-associated genes, indicating its critical role in mediating or regulating different genes that together contributes to the progression of DR. FCGR2B protein is a low affinity receptor for the Fc region of immunoglobulin gamma complexes.³⁸ FCGR2B involves many aspects of inflammatory and immune responses and the complex regulation of defense against infection.³⁹ However, the regulation of FCGR2B expression is complicated depending on the cell types. For example, interferon-gamma, a signature proinflammatory cytokine, escalates FCGR2B mRNA and cell surface expression by lipopolysaccharide-stimulated B cells while reducing it on monocytes, a process keeping immune response away from antibody production to pathogen clearance.³⁹ Increased levels of FCGR2B have been observed in the pancreata of autoantibody-positive at-risk individuals with type 1 diabetes versus controls,⁴⁰ and in the retina of a rat model with type 1 diabetes induced by streptocotocin.³¹ Similar to CCND1, FCGR2B is also one of the severity-associated genes with substantial interaction with other genes. Altogether, our results are in line with previous findings and further confirmed both CCND1 and FCGR2B increase in mRNA expression with the severity of DR, suggesting their potential role as biomarkers to monitor the progression of DR. 

In the early stages of DR, when blood-retinal barrier is intact, the immune system is mildly activated, and retinal microglia plays a major role in the para-inflammatory response, an important process to maintain retinal homeostasis.⁴¹ Consistent stimulation of damage-associated molecules induced by diabetes in the eyes cause maladaptation of the innate immune systems, ultimately promoting the development of DR.⁴¹ Circulating immune cells penetrate into the retina and result in chronic inflammation, leading to retina vascular and neuronal damage in the late stage.⁴² Despite the importance of the immune system in the pathogenesis of DR, few studies have investigated the changes in immune cell types over the progression of DR. Using the deconvolution method with CIBERSORTx, we systematically analyzed immune cellular profiles containing 22 immune cells and identified significant changes in proportions of resting mast cells, memory B cells, naïve B cells, and M2 macrophages in the macular region from donors with different degrees of DR. Mast cells are tissue-resident with heterogeneous phenotypes tuned by inflammatory stimuli.⁴³ We suspect that inflammation in DR results in the phenotype change of mast cells from resting to activated, for which proportions of resting mast cells are significantly decreased over the severity of DR. Memory B cells are the populations of cells that provide long-term humoral immunity.⁴⁴ A study revealed that chronic inflammation triggers an increase of activated memory B cells in patients with erythema nodosum leprosum, an inflammatory complication of leprosy.⁴⁵ Given that DR is a progressive disease featured with chronic inflammation, it is reasonable to see increased and decreased proportions of memory B cells and naïve B cells, respectively, with the severity of DR. Macrophages are heterogeneous and feature various functions once activated.⁴⁶ Our previous study demonstrated the significant increase of M2 macrophages in the retina of the rat model of retinal neovascularization and patients with PDR by deconvolution analysis.⁴⁷ Our results are in line with previous findings that the proportions of M2 macrophages are increased with the severity of DR. Correlation analyses showed CCND1 and FCGR2B have a moderate positive correlation with memory B cells and M2 macrophages, suggesting both genes may be crucial in regulating activation of M2 macrophages. 

Despite studies that have shown morphological or structural changes in retina cells from patients with DR, such as Müller glia and microglia cells,^48,49 few assessed changes in proportions of each retinal cell type from patients with DR. Using CIBERSORTx, we, for the first time, quantified the proportions of all retinal cells in the macular region from patients with different degrees of DR. We revealed that the proportion of Müller glia is significantly increased with the severity of DR, whereas changes of the proportions of other retinal cells were observed but without statistical significance. Expression levels of CCND1 and FCGR2B are moderate positive and related to the proportions of Müller glia, indicating their importance in the activation of Müller glia. Müller glia is important to maintain retinal function and health because it spans the entire width of the retina and contacts every other cell type. It is believed that Müller glia becomes activated in DR, characterized by the increased expression of glial fibrillary acidic protein, a marker of reactive gliosis.^49,50 However, some histological studies also reported Müller cells loss during the DR progress possibly due to the increased caspase-1 activity and IL-1β production following exposure to hyperglycemic conditions.^50,51 More studies are required to determine the mechanisms of Müller glia activation and death to pinpoint whether all Müller glia are equally affected by hyperglycemia. Decreased number of astrocytes in the retina has been found in both patients with DME and rat models of patients with diabetes.^52,53 Notably, our results showed an increased proportion of astrocytes, although it was not significantly correlated to the severity of DR. The DR progression is affected by multiple risk factors, such as age, gender, diabetes duration, fasting blood glucose, and glycosylated haemoglobin.^54,55 In our deconvolution analyses, we only controlled for age and gender as covariants of gene expression due to limited access to clinical labeling, leaving other risk factors that potentially confound the regression analyses. Moreover, a relatively small sample size could be another factor that affects the deconvolution outcome. Therefore, future studies with a large sample size and enriched clinical labeling should be performed to validate our results. 

There were several limitations of this study. First, the sample size in each group of the DR stage is relatively small, which likely affects the statistical power to identify a high number of severity-associated genes defined by LRT. In addition, the current study is based on in silico analysis of transcriptomic data. Future biological studies that investigate the functionality of CCND1 and FCGR2B in memory B cells and M2 macrophages, as well as Müller glia, would be important to understand their role in progressive DR. 

Conclusion

The findings reveal a group of genes, particularly CCND1 and FCGR2B, significantly increase in expression with the severity of DR, by controlling for other factors affecting gene expression, such as gender and age. Moreover, proportions of memory B cells, M2 macrophages, and Müller glia are significantly increased with the severity of DR. Our analyses provide potential avenues for subsequent research to study the mechanisms of DR progression. 

Code Availability

Analysis scripts will be made available upon request. 

Acknowledgments

The authors thank Ryan D. Chow for sharing codes for data analyses at Github. 

Supported by grants from the National Health and Medical Research Council of Australia (GNT1185600). The Centre for Eye Research Australia receives Operational Infrastructure Support from the Victorian Government. 

Author contributions: J.-H.W. and G.-S.L. conceived and designed the study. J.-H.W. developed the analysis approach, performed all data analyses, and generated figures. J.-H.W., R.C.B.W., and G.-S.L. prepared the manuscript. G.-S.L. supervised the work. 

Disclosure: J.-H. Wang, None; R.C.B. Wong, None; G.-S, Liu, None