Identification of potential molecular targets associated with proliferative diabetic retinopathy

Background This study aimed to identify and evaluate potential molecular targets associated with the development of proliferative diabetic retinopathy (DR). Methods The microarray dataset “GSE60436” generated from fibrovascular membranes (FVMs) associated with proliferative DR was downloaded from the Gene Expression Omnibus database. Differentially expressed genes (DEGs) from the active FVMs and control or inactive FVMs and control were evaluated and co-DEGs were identified using VEEN analysis. Functional enrichment analysis, and protein-protein interactions (PPI) network and module analyses were performed on the upregulated and downregulated coDEGs. Finally, several predictions regarding microRNAs (miRNAs) and transcription factors (TFs) were made to construct a putative TF-miRNA-target network. Results A total of 1475 co-DEGs were screened in active/inactive FVM samples, including 461 upregulated and 1014 downregulated genes, which were enriched for angiogenesis [Hypoxia Inducible Factor 1 Subunit Alpha (HIF1A) and Placental Growth Factor (PGF)] and visual perception, respectively. In the case of the upregulated co-DEGs, Kinesin Family Member 11 (KIF11), and BUB1 Mitotic Checkpoint Serine/Threonine Kinase (BUB1) exhibited the highest values in both the PPI network and module analyses, as well as the genes related to mitosis. In the case of downregulated co-DEGs, several G protein subunits, including G Protein Subunit Beta 3 (GNB3), exhibited the highest values in both the PPI network and module analyses. The genes identified in the module analysis were found to be from the signal transduction-related pathways. In addition, we were able to identify four miRNAs and five TFs, including miR-136 and miR-374. Conclusions In brief, HIF1A, PGF, KIF11, G protein subunits, and miR-136, miR-374 may all be involved in angiogenesis, retinal endothelial cell proliferation, and visual signal transduction in proliferative DR. This study provides a number of novel insights that may aid the development of future studies dedicated to discovering novel therapeutic targets in proliferative DR.


Background
Diabetic retinopathy (DR) is one of the most common microvascular complications in diabetes mellitus, and a primary cause of blindness and visual disturbance globally [1]. DR can be divided into non-proliferative and proliferative forms based on the type of microvascular lesions and related ischemic injury [2]. It has been predicted that the global prevalence of diabetes will continue to grow significantly for the next several decades, as a result of increasing incidence of type 1 diabetes rather than increasing type 2 diabetes which has already reached an epidemic level [3]. The incidence of DR increases with diabetic progression and results from a complex process with various molecular and biochemical participants, including oxidative stress, endoplasmic reticulum stress and mitochondrial damage amongst others [4]. The basic mechanisms underlying DR have been extensively studied, however, effective prevention and therapeutic interventions for this disease are still not available.
At present, DR treatments mainly focus on laser photocoagulation, intravitreal pharmacotherapy, anti-VEGF, and glycemic control [5]. Intravitreal anti-inflammatory agents have been effective in treating DR [6], and intravitreal triamcinolone acetonide shows an antiangiogenic effect toward proliferative DR [7]. Vascular endothelial growth factor (VEGF) plays a critical role in promoting vascular permeability, cell migration as well as proliferation of vascular endothelial cells, vasculogenesis, and angiogenesis [8]. VEGF is intimately involved with the progression of proliferative DR and diabetic macular edema facilitating changes in retinal capillary permeability, and advances in anti-VEGF therapies, in age-related macular degeneration, have accelerated the application of anti-VEGF therapies in DR [9]. Oxidative stress is a critical factor in the etiology of DR. Metabolic abnormalities, caused by the increase of glucose concentration in diabetes, can lead to the overproduction of the superoxide radical involved in the uncoupling of mitochondrial electron transport chains, ultimately resulting in oxidative stress [10,11]. Antioxidants, including Vitamins and polyphenols, are considered a beneficial therapeutic strategy during the treatment of DR. These compounds inhibit the reactive oxygen species, free radicals and enhance the antioxidant defense system [12]. Despite some progress, DR remains a prevalent vision-threatening disease.
In this study, we used the microarray dataset "GSE60436" generated by Ishikawa K et al., [13]. In that study, they identified a subset of differentially expressed genes (DEGs) from active fibrovascular membranes (FVMs) and normal retinas or inactive FVMs and normal retinas, these DEGs were then evaluated by functional analysis. In our study, we aim to identify pivotal genes associated with the pathogenesis of proliferative DR and explore their function and upstream and downstream targets. The common DEGs from active/ inactive FVMs and normal retinas were identified, and then subjected to functional enrichment, proteinprotein interaction (PPI) network and module analyses. In addition, the miRNA-target and transcription factors (TFs)target were predicted to explore the potential regulatory relationships.

Study approval
This study did not use any animal or human participants. All data was generated from public databases.
Limma package software (Version 3.10.3, http://www. bioconductor.org/packages/2.9/bioc/html/limma.html) was used to process the raw CEL files downloaded from GEO, and data preprocessing was done using the RMA (robust multi-array average) method including background correction, normalization, and expression calculations. The probes were removed when they were not able to be matched to a specific gene symbol, and the average value was taken as the expression value for each gene when different probes matched to the same gene symbol.

DEG screening
The genes that were differentially expressed in different groups [active-FVMs vs. control (active group); inactive-FVMs vs. control (inactive group)] were analyzed using the empirical Bayes test from the limma package software [14]. The DEGs were screened using a cut-off value of P < 0.05 and |log fold change (FC)| > 2. The DEGs from both groups were then subjected to VEEN analysis using the VENNY (Version 2.1.0, http://bioinfogp.cnb. csic.es/tools/venny/index.html) online tool. Overlapping genes were considered co-regulated DEGs (co-DEGs) in the following analyses.

Functional enrichment analysis
The Database for Annotation, Visualization and Integrated Discovery (DAVID, Version 6.8, https://david-d. ncifcrf.gov/) online tool was used to analyze the KEGG pathway and Gene Ontology annotations for the co-DEGs. The number of enriched genes was set as: count ≥2, and P < 0.05 was considered the threshold value for significantly enriched terms.

PPI network and module analysis
Search Tool for the Retrieval of Interacting Genes (STRING, Version 11.0, http://www.string-db.org/) was used to predicted the interactions among the co-DEGs with the parameters set to species = human, and PPI score = 0.9 (highest confidence). The PPI network was constructed using Cytoscape software (version 3.2.0, http://www.cytoscape.org/) based on the interactions retrieved from STRING. In addition, the MCODE plugin [15] (Version 1.4.2, http:// apps.cytoscape.org/apps/MCODE) from Cytoscape was used to identify the significant modules, with scores > 10, and was followed by functional enrichment analysis.

TF-miRNA-target network construction
The overrepresentation enrichment analysis method from WebGestalt [16] (http://www.webgestalt.org/) was used to predict the TF-and miRNA-target interactions for all of the genes from the significant modules using a cut-off value of P < 0.05. Then, the TF-miRNA-Target regulatory network was constructed by combining miRNA-target and TF-target regulatory interactions.

Data preprocessing and DEG screening
The expression values for all the genes from the nine samples were normalized using the RMA method, and values with an unchanged position in the boxplot were used for subsequent analysis, as this can be used as a proxy for normalization (Fig. 1a). A total of 2025 DEGs were identified in the active group, including 758 upregulated DEGs and 1297 downregulated DEGs. Similarly, 1961 DEGs (757 upregulated and 1204 downregulated) were identified in the inactive group. Fig. 1b shows the heat maps of those DEGs, and reveals that the DEGs can be easily distinguished from each of the samples. In addition, the VEEN analysis, which was performed to identify common DEGs in the two groups, identified 1475 overlapping or co-DEGs. Of these 461 genes were upregulated, and 1014 genes were downregulated (Fig. 1c).

Functional enrichment analysis
The upregulated co-DEGs were significantly enriched for 148 GO-biological processes (GO-BP) and 41 KEGG pathways. Including "GO:0030198~extracellular matrix organization", "GO:0001525~angiogenesis" and pathways "hsa04512: ECM-receptor interaction", "hsa04151: PI3K-Akt signaling pathway". In addition, the downregulated co-DEGs were significantly enriched for 150 GO-BP Fig. 1 Data normalization and the distribution of differentially expressed genes. a, Box plots illustrating data normalization. b, Heat maps of DEGs in active and inactive groups. The horizontal axis represents each sample, and the left vertical axis shows clusters of DEGs. The gradual change of color from green to red represents changes in the expression values from low to high. c, Venn diagram of DEGs in the two groups. DEGs, differentially expressed genes terms and 40 KEGG pathways including, "GO: 0007601~visual perception", "GO:0001523~retinoid metabolic process" and pathways "hsa04744: Phototransduction" and "hsa04724: Glutamatergic synapse". Table 1 summarizes the top 10 GO-BP terms and KEGG pathways.

PPI network and module analysis
For the upregulated co-DEGs.
The upregulated co-DEGs were used to retrieve various predicted interactions, which were then used to create the PPI network. This network consisted of 202 nodes and 751 interactions (Supplemental Fig. 1). Moreover, two significant modules with score ≥ 10 were identified. Module-A (score = 18.4) contained 21 nodes and 184 interactions, and module-B (score = 12) included 12 nodes and 66 interactions (Fig. 2). Kinesin Family Member 11 (KIF11), BUB1 Mitotic Checkpoint Serine/Threonine Kinase (BUB1) and Cyclin B2 (CCNB2) were the hub genes with the highest degree of relevance in both the PPI network and module-A. The genes in module-A were found to be related to three KEGG pathways and 25 GO-BP terms, including "GO:0007067~mitotic nuclear division" and "GO: 0051301~cell division". Additionally, the genes in module-B were found to participated in three KEGG pathways and eight GO-BP terms, including "GO:0030198~extracellular matrix organization".
For the downregulated co-DEGs. The PPI network constructed using downregulated co-DEGs contained 349 nodes and 870 interactions (Supplemental Fig. 2), from which two significant modules with score ≥ 10 were also identified. Module-A (score = 18.24) contained 26 nodes and 228 interactions and module-B (score = 10) included 10 nodes and 45 interactions (Fig. 3). G protein subunits were the hub genes with the highest degree of conservation in both the PPI network and module-A, including G Protein Subunit Beta 3 (GNB3), GNB5, G Protein Subunit Gamma 13 (GNG13), GNG8, GNG17, and GNG3. The genes in module-A were found to be significant in 16 KEGG pathways and 34 GO-BP terms, including "hsa04725: Cholinergic synapse" and "GO:0007186~G −protein coupled receptor signaling pathway". Similarly, the genes in module-B were found to be associated with two KEGG pathways and nine GO-BP terms. Fig. 4 shows the enriched KEGG pathways and top five GO-BP terms for the genes found in the significant modules.

TF-miRNA-target network construction
The genes from the significant modules were used to predicted the TF-and miRNA-target interactions, which were then used to construct the TF-miRNA-target network (Fig. 5). The TF-miRNA-target network consisted of 33 nodes and 27 interactions, from which we were able to identify four miRNAs and five TFs, which may interact with up to 24 novel genes. GNG3 was shown to be regulated by miR-136, and Growth Factor Independent 1 Transcriptional Repressor (GFI1). KIF20A was regulated by miR-374, and so on.

Discussion
In this study, we identified 1475 co-DEGs from active/inactive-FVMs samples. The upregulated co-DEGs were found to participate in angiogenesis [Hypoxia Inducible Factor 1 Subunit Alpha (HIF1A) and Placental Growth Factor (PGF)] and mitosis-related processes (the genes in upregulated module-A, including KIF11 and BUB1). The downregulated co-DEGs were enriched in G − protein coupled receptor signaling and other pathways (The genes in downregulated module-A, including GNB3 and GNG8).
The important pathological features of proliferative DR including abnormal growth of retinal blood vessels and angiogenesis have all been linked to VEGF signaling. High blood sugar can trigger hypoxia in retinal tissues, and this hypoxia acts as a crucial factor in the regulation of VEGF-induced angiogenesis through the production of HIF [17]. HIF1A is an hypoxia induced TF and the lack of HIF1A suppresses the formation of the retinal intermediate vascular plexus in mouse models [18]. PGF, a homologue of VEGF, plays an important role in placental development, and has been shown to be involved in pathological angiogenesis both in ocular and nonocular cancers [19]. In addition, increased expression of PGF has also been observed in the vitreous of DR patients [19]. Bender et al., revealed that PGF and VEGFA are both essential for follicular angiogenesis in primates, and neither can function alone and result in normal follicular angiogenesis [20]. Huang et al. showed that PGF negatively effects retinal endothelial cell barrier function by activating VEGFR1 and VEGFR2 and inhibiting glucose-6-phosphate dehydrogenase and the antioxidant pathways [21]. Ishikawa K et al., reported that the genes that were significantly upregulated in active FVMs were predominantly part of the angiogenesis process [13]. Therefore, we speculate that PGF and VEGF are important in the progression of FVM to proliferative DR.
MicroRNAs (miRNAs) are small non-coding RNAs of around 22 nucleotides (nt) in length, which mediate the expression of protein-coding genes, and are becoming increasingly popular biomarkers used in both prognostic and diagnostic assays [22,23]. Targets for miR-136 and miR-374 were significantly enriched in the co-DEGs identified in this study. However, the relationship between miR-136 or miR-374 and DR is not well reported, and this is the first publication linking them. miR-136 is highly expressed in pre-eclampsia and inhibits the formation of HUVEC capillaries via dysregulation of VEGF. In addition, the expression levels of miR-136 inversely correlates with VEGF [24]. In paclitaxel-resistant ovarian cancer cells, pre-miR-136 transfection significantly decreased angiogenesis and induced   apoptosis, and regulation of miR-136 expression could resensitize paclitaxel-resistant cells [25]. C-C Motif Chemokine Ligand 3 represses the expression of miR-374b, which subsequently accelerates VEGF-A level and angiogenesis in osteosarcoma cells [26]. Thus, we speculate that miR-136 and miR-374 may participate in the regulation of angiogenesis in proliferative DR. KIF11 is a member of the kinesin-like protein family. Mutations in KIF11 have been shown to cause familial exudative vitreoretinopathy associated with retinal detachment [27]. Birtel et al., revealed that KIF11 is important in both ocular development and the maintenance of retinal morphology and function and defects in this protein have been associated with retinal ciliopathy [28]. KIFs are motor proteins involved in the activities of the centrosome and spindle during mitosis. In our study, KIFs (KIF11, KIF15, KIF20A) and other mitosis-associated genes, like BUB1, were significantly enriched in one module involved in cell division and mitotic nuclear division. The proliferation of retinal endothelial cells is a significant event in the progression of DR [29]. In addition, a previous study has shown that high glucose significantly increases the proliferation of retinal endothelial cells and the expression levels of VEGF, also, it was more effective under intermittent high glucose, than constant high glucose conditions [30]. Therefore, we suggest that the genes in this module play a key role in proliferative DR through the regulation of the proliferation of retinal endothelial cells. Fig. 2 The significant modules identified in the PPI network of upregulated co-DEGs. Node size represents the degree score; lines represent interactions; co-DEGs, co-regulated differentially expressed genes Fig. 3 The significant modules identified in the PPI network of downregulated co-DEGs. Node size represents the degree score; lines represent interactions. co-DEGs, co-regulated differentially expressed genes  G proteins are a kind of heterotrimeric guanine nucleotide-binding protein consisting of three subunits (alpha, beta and gamma subunits), which play an important role in intracellular signal transduction, including insulin, epinephrine, dopamine signaling molecule transduction [31]. GNB3 encodes G Protein Subunit Beta 3. The GNB3 gene polymorphisms, including C825T, have been linked to type 2 diabetes mellitus and its complications [32]. GNB3 knockout lead to the malfunction of cone photoreceptors and ON-bipolar cells in murine retinas [33]. In addition, it was revealed that a naturally occurring mutation in GNB3 results in retinal degeneration in chickens [34]. A homozygous missense variant in GNB3 results in a unique stationary retinal disorder with dual anomalies in visual processing, night blindness and photophobia [35]. G protein subunits (GNB3, GNB5, GNG3, GNG7, GNG8, GNG13) are downregulated in proliferative DR, and significantly enriched in module A which is associated with cellular responses to glucagon stimuli including G − protein coupled receptor signaling pathways and other biological processes. Therefore, we conclude that G proteins may play a crucial role in visual signal transduction in proliferative DR.
Although this study has made several novel observations, all the results were based on bioinformatics analysis, and further experimental verification is required. In addition, in an ideal situation, healthy retinas, nonproliferative diabetic retinas and autologous retinas should be used as controls for a more comprehensive comparison with retinas from subjects with FVM-associated proliferative DR. However, such microarray studies are limited by challenges in obtaining samples. Additionally, the small sample size in this study may also limit the application of our observations. Genes that are important for proliferative DR might be excluded as a result of the small sample size and inter-sample variation.

Conclusions
In brief, we investigated differentially expressed gene profiles in proliferative DR. Several genes, including HIF1A, PGF, KIF11, G protein subunits, and miR-136, miR-374 were found to be involved in angiogenesis, retinal endothelial cell proliferation, and visual signal transduction in proliferative DR. These findings could provide valuable insight and foundations for the study of proliferative DR's pathomechanism and may provide an idea of potential therapeutic targets for its treatment.
Additional file 1 Figure S1 The PPI network of upregulated co-DEGs. Node size represents the degree score; lines represent interactions; co-DEGs, co-regulated differentially expressed genes.
Additional file 2 Figure S2. The PPI network of downregulated co-DEGs. Node size represents the degree score; lines represent interactions; co-DEGs, co-regulated differentially expressed genes.