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.

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.

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.

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.

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.

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.