Hub Genes Associated with the Diagnosis of Diabetic Retinopathy

Introduction

Diabetes mellitus (DM) affects an increasing number of people worldwide.¹ The global incidence of DM is estimated to be 415 million adults and is expected to be nearly 700 million by 2045.^2,3 The course of DM is closely related to the development and progression of diabetic retinopathy (DR),⁴ which is one of the most common complications of DM and which can reduce visual acuity and even lead to blindness.^5,6 Almost all patients with type 1 diabetes mellitus (T1DM) and more than 60% with type 2 diabetes mellitus (T2DM) have some degree of retinopathy after 20 years with DM.⁷ The estimated global prevalence of DR is 382 million,⁸ imposing a heavy financial burden on healthcare systems.^9,10

The retina is the most active metabolic tissue and consumes a large amount of energy and oxygen.¹¹ Inadequate oxygen and nutrient supply to the retina may impair the oxygen balance in cells, resulting in a hypoxic retinal environment.¹² During the development of DM, high blood glucose reduces the oxygen levels in the retina,¹³ while chronic hyperglycemia alters the function of vascular smooth muscle cells, thus promoting DR.¹⁴ A typical manifestation of DR is gradual change in retinal microvasculature, which in turn leads to retinal ischemia, neovascularization, changes in retinal permeability, and macular edema.¹⁵ Despite the existing knowledge about the pathology and manifestation of DR, no prophylactic treatment has been reported for the early stages of the disease.¹⁶

DR can be clinically divided into non-proliferative diabetic retinopathy (NPDR), which is characterized by microvascular damage, and proliferative diabetic retinopathy (PDR), where new blood vessels are formed.¹⁷ NPDR can also develop into PDR, resulting in irreversible vision loss and complete blindness.¹⁸ Current treatments for DR include laser surgery and anti-vascular endothelial growth factor (VEGF) therapy, which can control the growth of abnormal blood vessels but cannot prevent blindness.¹⁹ Gene therapy can also be used as a restorative or neuroprotective approach depending on when it is applied during the disease and what genes it targets.²⁰ DNA methylation in DM patients has also been associated with the pathogenesis of DR, suggesting that methylated sites can serve as predictive or diagnostic biomarkers of DR.²¹

Here, we used the GSE146615 dataset to identify potential diagnostic markers and therapeutic targets in DR, as well as to begin to clarify molecular dysregulation in the disease. In addition, we used the GSE76169 dataset to explore potential associations between methylation mechanisms and DR. In GSE76169 dataset, all T1DM cases received conventional treatment in the Diabetes Control and Complications Trial (DCCT)²² and developed retinopathy at the 10th year after the end of the EDIC. In contrast, controls were in the DCCT intensive treatment group and did not develop retinopathy during the EDIC study. Identifying early features of DR is particularly important for inhibiting pathology and preventing visual impairment. Therefore, we expect that understanding the underlying pathogenic mechanisms of DR will contribute to the development of novel therapeutic approaches.

Materials and Methods

Data Processing

The DR gene data were collected from the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo). The GSE146615 dataset included gene expression profiles of lymphoblastoid cell lines derived from healthy controls (NoD) and matched T1DM individuals with retinopathy (DwC) and without retinopathy (DwoC). The cells were cultured under standard glucose (SG) (21 NoD, 26 DwC, and 25 DwoC) and high glucose (HG) conditions (21 NoD, 26 DwC, and 25 DwoC), and all samples were analyzed on Illumina HT12v4 microarrays (ThermoFisher).

The raw data were background-subtracted and normalized using the lumiR package.²³ The GSE76169 dataset included the whole-blood DNA methylation profiles from 32 DR cases and 31 controls collected at the beginning of the Epidemiology of Diabetes Interventions and Complications (EDIC) study.²⁴

DNA methylation samples were measured on the Illumina HumanMethylation450 BeadChip array, and the Bioconductor package ChAMP was used to generate the normalized beta values.²⁵ All data in the present study were obtained from an open-access database and no ethical approval was required.

Differential Analysis

The limma package inR²⁶ was used for the differential analysis between DR and DM patients or between DM patients and controls in the GSE146615 dataset under HG or SG conditions. Genes with P value < 0.05 of t test were considered statistically significant and were indicated as differentially expressed genes (DEGs). In contrast, common genes were defined as genes whose expression in DR differed in the same direction relative to DM under both HG and SG conditions. The limma algorithm was also applied to the GSE76169 dataset to obtain the differentially methylated CpG sites with P < 0.05.

Network Construction

The Search Tool for Retrieval of Interacting Genes/Proteins (STRING; http://string-db.org) is an online database for analyzing functional networks between two or more proteins. Here, the protein–protein interaction (PPI) network was constructed by analyzing common genes with STRING based on the criterion of a combined score >500. The network was displayed using Cytoscape software. Significant modules in the PPI network were identified by the molecular complex detection (MCODE) plugin in Cytoscape²⁷ with degree cutoff = 2, node score cutoff = 0.2, K-Core = 5, and max. depth = 100.

Enrichment Analysis

Annotation of signaling pathways and functions for module genes were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID, http://david.ncifcrf.gov). Genes were also analyzed for enrichment of Gene Ontology (GO) processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways pathways using the clusterProfiler package in R.²⁸ The same package was also used within a Java environment to obtain the background set for gene set enrichment analysis (GSEA) between DR and DM.²⁹ The fast preranked gene set enrichment analysis (fgsea) package inR³⁰ was used to further analyze the GSEA results. Differences with P value < 0.05 were considered statistically significant. Gene set variation analysis (GSVA) was carried out using the GSVA package³¹ to explore the activation or inhibition of biological processes and KEGG pathways in DR relative to DM.

Short Time-Series Expression Miner Analysis

Short Time-series Expression Miner (STEM) software³² was used to cluster module genes of the GSE146615 dataset in healthy, DwoC, and DwC samples under SG and HG conditions. A P value < 0.05 was considered as statistically significant clustering. The expression of the significantly clustered genes was progressively up- or down-regulated in the sequence: healthy < DwoC < DwC. The area under the receiver operating characteristic curve (AUC) was also calculated to evaluate the diagnostic value of the hub genes using the pROC package.³³

Single-Sample GSEA

The marker gene sets defining different immune cell types were obtained from a previous study.³⁴ Differences in infiltration levels by immune cells between DR and DM was calculated using the single-sample gene set enrichment analysis (ssGSEA) function in the GSVA package and was expressed as a ssGSEA score. A P value < 0.05 was considered statistically significant.

Results

Identification of DEGs Between DR and DM

In order to identify specific genes associated with DR, differential analysis of gene expression between DR and DM in GSE146615 was first performed (Figure 1), resulting in 3155 DEGs under HG conditions (Figure 2A, Table S1) and 3045 DEGs under SG conditions (Figure 2B, Table S2). The intersection of DEGs were attributed to the progression of DR and indicated that 1053 common genes were upregulated and 864 were downregulated under the two glucose conditions (Figure 2C).

Figure 1 Study flowchart. The GSE146615 dataset was used to identify potential diagnostic biomarkers and therapeutic targets of diabetic retinopathy (DR), as well as molecular mechanisms of the disease. Alterations in methylation of DR-related genes were also screened to identify methylated markers. Abbreviations: HG, high glucose; SG, standard glucose; PPI, protein–protein interaction; STEM, Short Time-series Expression Miner; MCODE, molecular complex detection; GSVA, gene set variation analysis; AUC, area under the receiver operating characteristic curve.

Figure 2 The differentially expressed genes in DR compared to DM. Differentially expressed genes between diabetic retinopathy (DR) and diabetes mellitus (DM) under (A) high glucose (HG) conditions and (B) standard glucose (SG) conditions. Red indicates upregulated genes and green indicates downregulated genes. The top eight genes with the largest fold change are listed. (C) A total of 1053 differentially expressed genes that were upregulated under both HG and SG conditions were detected, as well as 864 genes that were downregulated under the two glucose conditions. These DEGs were defined as “common genes”.

PPI Network

The protein interactions among the identified common genes were predicted using the STRING database. In particular, the 1917 common genes were searched in the STRING (http://string-db.org) and screening for interaction scores >500 yielded a network of 6775 reciprocal pairs of PPIs. In addition, Cytoscape software was used to identify the interaction relationship of the PPI network genes, and the MCODE plugin recognized 279 module genes (Figure 3A). Further calculation of their AUC values under HG and SG conditions indicated that the following eight genes had AUC > 0.8: FOXO1, ALDH9A1, TSC1, AHI1, CMKLR1, EXOSC7, INTS6, and MDFIC. These eight genes were identified as hub genes (Figure 3B). The differential expression of these genes between DR and DM was similar under SG conditions (Figure 3C) and HG conditions (Figure S1).

Figure 3 Construction of the protein–protein interaction (PPI) network and identification of hub genes. (A) PPI network of eight modules. Differently colored nodes represent different modules. (B) AUC values of the eight module genes under high glucose (HG) and standard glucose (SG) conditions (left). Mean AUC for genes with AUC > 0.8 under HG and SG conditions (right). Red and blue nodes represent, respectively, genes upregulated or downregulated in diabetic retinopathy (DR) relative to diabetes mellitus (DM). (C) Differential expression profiles of hub genes between DR and DM under SG conditions. ***P < 0.05. (D) Differences in immune cell infiltration between DR and DM under SG conditions. Each row represents a type of immune cell. Blue lines represent significant downregulation in DR. (E) Correlation between hub genes and infiltrating immune cells based on Pearson correlation analysis under SG conditions. Red blocks represent positive correlations; blue, negative correlations. *P < 0.05, **P < 0.01. Abbreviation: AUC, area under the receiver operating characteristic curve.

The infiltration by immune cells in DR and DM was also quantified, revealing differences between the two conditions. Th2 and B cells were significantly downregulated in DR compared to DM under both SG (Figure 3D) and HG (Figure S2) conditions. Correlation analysis also showed that infiltration by Th2 cells correlated positively with AHI1 expression and negatively with TSC1 expression under both glucose conditions (Figure 3E and Figure S3). The infiltration by B cells correlated positively with MDFIC expression and negatively with FOXO1 expression.

Enrichment Analysis of Module Genes

The module genes were studied by enrichment analysis to identify the molecular mechanisms of DR. There were 414 terms in the DAVID analysis results, included 5′‐adenosine monophosphate‐activated protein kinase (AMPK) signaling pathway, cell cycle, and tumor necrosis factor (TNF) signaling pathway (Table S3). In the enrichment results, the activation or inhibition of biological processes and KEGG terms in DR relative to DM was examined under SG conditions using the GSVA package. Specifically, vascular smooth muscle cell development, regulation of retinoic acid biosynthetic process, and regulation of oxygen metabolism were activated in DR and significantly enriched in module genes (Figure 4A). In contrast, the activation of c-Jun N-terminal (JUN) kinase, mitogen‑activated protein kinase (MAPK), and chronic inflammatory response were inhibited in DR. The KEGG pathways of valine, leucine, and isoleucine degradation, AMPK signaling, and β-alanine metabolism were also activated in DR and significantly enriched in module genes (Figure 4B), whereas the cell cycle and signaling pathways of Toll-like receptor (TLR) and TNF were inhibited. Furthermore, GSEA indicated the activation of histidine metabolism, maturity-onset diabetes of the young, and neutrophil extracellular trap formation under SG (Figure 4C) and HG (Figure 4D) conditions in DR. GSEA also indicated inhibition of the T-cell receptor, TLR, and NFκB signaling pathways in DR.

Figure 4 Module genes involved in biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. (A and B) Numbers of module genes in diabetic retinopathy compared to healthy controls involved in activation or inhibition of (A) biological processes and (B) KEGG pathways, as determined by GSVA. Longer columns indicate the involvement of more genes under the given term. (C and D) Module genes involved in activation or inhibition of KEGG pathways in diabetic retinopathy compared to healthy controls under (C) standard glucose or (D) high glucose conditions, as determined by GSEA. Differences with P < 0.05 were considered statistically significant. Abbreviations: GSVA, gene set variation analysis; GSEA, gene set enrichment analysis; BP, biological processes; NES, normalized enrichment score.

Methylation of Common Genes

The methylation profiles of DR cases and controls in the GSE76169 dataset were studied by differential analysis to determine if the expression of common genes might be regulated by methylation. Genes were considered to be methylation markers regulated by methylation when their expression differed (under SG conditions) in the opposite direction to methylation differed in DR patients compared controls. A total of 2992 differentially methylated probes (DMPs) were identified (Figure 5A) and after comparing them with common genes, we found 92 methylation markers that may be involved in the progression of DR (Figure 5B).

Figure 5 Identification of methylation markers for common genes. (A) Differentially methylated positions between diabetic retinopathy and diabetes mellitus. Among them, 71.19% were hypermethylated and 28.81% were hypomethylated. Different colors represent different methylation sites. Labels indicate the 10 sites with the highest methylation levels. (B) Methylation and expression levels of methylation markers. The labels indicate the important genes. (C) Biological functions and (D) KEGG pathways significantly enriched by methylation markers. Thicker lines represent greater correlation between the terms and methylated markers. Abbreviations: DMPs, differentially methylated positions; KEGG, Kyoto Encyclopedia of Genes and Genomes; corr, Pearson correlation coefficient.

Enrichment analysis revealed that these methylation markers were significantly involved in biological processes, such as “vascular smooth muscle cell development”, “reactive oxygen species biosynthetic process” (Figure 5C), In KEGG results, we identified significantly enriched “Notch signaling pathway”, “Glycine”, “serine and threonine metabolism” (Figure 5D).

Potential Genes Involved in DR Progression

Through STEM analysis, we identified the genes whose persistent expression changes during the progression from healthy controls to DM patients and then to DR patients, under both HG and SG conditions. A total of 18 and 17 common genes showed progressively up- or down-regulated expression under SG and HG conditions, respectively, and the groups of genes fell into two clusters (Figure 6A and B). Among the two groups of progressively dysregulated genes, three intersection genes were found to be module genes of the PPI network (CXCL10, HCK, and HLA-DRB1), and were therefore considered as key genes. The expression of key genes was progressively downregulated from healthy controls to DM patients and then to DR patients, under both SG (Figure 6C) and HG (Figure 6D) conditions.

Figure 6 Identification of progressively dysregulated genes during the development of diabetic retinopathy. (A and B) Progressively dysregulated genes under (A) standard glucose (SG) and (B) high glucose (HG) conditions were identified by STEM analysis in the trend: healthy controls < diabetes mellitus < diabetic retinopathy. The genes were then clustered into a progressively upregulated cluster (Cluster 1) and a progressively downregulated cluster (Cluster 2). (C and D) Expression of CXCL10, HCK, and HLA-DRB1 genes among healthy controls, diabetic retinopathy, and diabetes mellitus under (C) SG and (D) HG conditions. Abbreviations: STEM, Short Time-series Expression Miner; NoD, healthy controls; DwC, type 1 diabetes individuals with retinopathy; DwoC, type 1 diabetes individuals without retinopathy.

Discussion

The identification of effective biomarkers for DR is particularly important, as such biomarkers may contribute to the early diagnosis and treatment of the disease.¹⁵ Although biomarkers with therapeutic potential are commonly used in medicine and especially in oncology,^35,36 none has been reported for DR. In the present study, we investigated potential diagnostic and therapeutic biomarkers of DR by comparing gene expression changes in DR and DM cell lines under SG and HG conditions.

Screening of genes with different expression in the same direction under SG and HG conditions revealed that several genes may play a key role in DR diagnosis. In contrast to DM, FOXO1, ALDH9A1, and TSC1 are upregulated in DR, thus promoting the disease. Specifically, the upregulated expression of FOXO1 causes IL-1β-induced autostimulation,³⁷ and its upregulation in endothelial cells significantly contributes to the development of DR.³⁸ In addition, FOXO1 acts downstream of insulin-dependent signaling pathways, which are dysregulated in DM.³⁹ TSC1 is another important gene that contributes to the development of diabetic complications, including DR, by targeting the mechanistic target of rapamycin (mTOR) pathway.^40,41 Among the genes downregulated in DR, AHI1 encodes a cilia-localized protein, and its absence promotes retinal degeneration.⁴² CMKLR1 is activated by the Wnt/β-catenin pathway, which is involved in extracellular matrix accumulation and renal fibrosis in DR.^43,44 However, the pathogenic role of CMKLR1 in retinal microvascular function has not been evaluated in vivo or in vitro. EXOSC7 has been reported to be differentially expressed between T1DM and healthy controls,⁴⁵ but no study has been found showing its altered expression in DR. INTS6 is a glucocorticoid response gene, and glucocorticoids have gluconeogenesis promoting effects on the liver.⁴⁶ Levels of INTS6 are significantly upregulated after intermittent fasting in patients with metabolic syndrome.⁴⁷ The glucocorticoid receptor binds to MDFIC, and directly inhibits its expression.⁴⁸ And MDFIC was significantly downregulated in retinopathy compared with the control.⁴⁹ Although it is clear that the up- and down-regulation of hub genes plays a significant role in the DR process, their specific roles in the disease remain to be elucidated. In addition, the hub genes may affect the immune response of DR, according to our correlation results between immune cell infiltration and gene expression.

Here, we first performed enrichment analysis of module genes within the PPI network and identified potential biological functions and signaling pathways that may be involved in the pathophysiology of DR. Oxidative metabolic processes are strongly activated in DR cells, as hypoxia in the retina can induce neuronal degeneration.⁵⁰ Alterations in the excitability of vascular smooth muscle cells play a key role in the development and progression of vascular complications in DM.^51,52 Animal studies have shown that hyperglycemia can cause vascular dysfunction in DR,⁵³ while inflammation has also been indicated as a key factor in the development of DR.⁵⁴ However, here, we found that chronic inflammation is suppressed during DR, suggesting that the inflammatory response of in vitro cultured DR cells may not play a major pathological role in the disease development.

Enrichment analysis of KEGG terms showed that the AMPK signaling pathway is activated in DR. This pathway mediates the stress response, promotes autophagy, and stabilizes hypoxia-inducible factor 1α (HIF-1α), thus enhancing VEGF expression.⁵⁵ VEGF hypersecretion then leads to irregular retinal vasculature, resulting in edema, hemorrhage, and eventually retinal detachment.^56,57 However, we found that inflammatory signaling pathways, such as those involving TLR and TNF, are suppressed in DR, which is inconsistent with previous studies that suggested that HG levels enhance the expression of TLR, NFκB, and TNF-α, thereby exacerbating DR.⁵⁸ Our results suggest that DR cells cultured in vitro may be more prone to activate oxidative stress-related signals, while suppressing inflammation-related responses. This may be related to the reduced immune cell content in DwC compared to DwoC of GSE146615 dataset.

STEM analysis showed that the expression of CXCL10, a known biomarker for DR,⁵⁹ as well as expression of HCK and HLA-DRB1 was progressively downregulated from healthy controls to DM patients and then to DR patients. This suggests that the three genes may be involved in the progression of DR in patients with DM. These results are inconsistent with recent studies, where CXCL10 expression was upregulated from DM to DR.⁵ This may be related to the suppression of inflammatory response in DR in our results. Platelet-derived vesicles seem to regulate hyperglycemia-induced retinal endothelial damage by releasing CXCL10, which can then activate the TLR4 pathway.⁶⁰ HCK belongs to the Src family kinase involved in regulating endothelial cell proliferation, migration, and angiogenesis,⁶¹ and it is differentially expressed between DM and normal subjects, as well as between healthy controls and patients with diabetic peripheral neuropathy.^62,63 However, its role in DR progression has not yet been elucidated. HLA-DRB1, for its part, is a vision-related gene that is differentially expressed in proliferative DR lesions,⁶⁴ while its polymorphisms are known risk factors for DM and its complications.^65,66 Based on these earlier findings and our knowledge of DR pathophysiology, we suggest that CXCL10, HCK, and HLA-DRB1 play an important role in the development and progression of the disease.

Differential analysis of the methylation profiles of DR cases and controls in the GSE76169 dataset identified altered DNA methylation of DR-related genes. Consistent with our results, transcription of genes related to redox balance is known to be regulated by cytosine methylation status in DNA. This methylation is involved in retinal development and responds dynamically to environmental changes.⁶⁷ Indeed, the restoration of proper DNA methylation can restore mitochondrial homeostasis, maintain cellular redox balance, and inhibit or delay DR progression.⁶⁸ In the KEGG pathways of methylation markers enriched, the balance of Notch signaling controls the formation and degeneration of retinal vessels.⁹ H₂O₂ caused a decrease in the concentration of glycine in retina, which is involved in the process of DR.⁶⁹ Hyperglycemia induced activation of serine/threonine kinase family proteins is one of the factors involved in microvascular alterations in DM that can lead to DR.⁷⁰

Our results should be interpreted with caution. The small sample in our study means that our findings should be verified and extended in larger studies. Moreover, our findings were not validated through molecular or clinical experiments, which will be performed in subsequent studies to confirm the efficacy of the identified genes as novel diagnostic biomarkers or therapeutic targets in DR. In addition, although T1DM and T2DM overlap in some clinical and pathological features, their molecular mechanisms and genetic basis are different. We only selected T1DM individuals to do the analysis, whether these results are applicable in T2DM needs further study.

Conclusion

Our study on differentially expressed genes between DR and DM showed that oxidative stress may be an important pathological mechanism in DR cells cultured in vitro. A total of eight genes were identified that could be used as valuable biomarkers for DR, while CXCL10, HCK, and HLA-DRB1 could serve as biomarkers of DR severity.