Multiple Bioinformatics Analyses of Integrated Gene Expression Profiling Data and Verification of Hub Genes Associated with Diabetic Retinopathy.

BACKGROUND Diabetic retinopathy (DR) is a serious complication of diabetes that can lead to blindness. This study aimed to identify the core genes and molecular functions involved in DR through multiple bioinformatics analyses. MATERIAL AND METHODS The mRNA gene profiles of human DR tissues from the GSE60436 and GSE53257 datasets were assessed with R software and integrated to identify the co-expressed differentially expressed genes (DEGs). Multiple bioinformatics analyses were used: Gene Ontology (GO) analysis, signaling pathway analysis, and hub gene prediction. Quantitative reverse transcription-PCR (qRT-PCR) was used to verify the hub genes. RESULTS The Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool suggested that the biological processes of the DEGs focused on mitochondrial transport, the cellular components focused on mitochondria, and molecular functions focused on catalytic activity. The results provided by DAVID were consistent with those provided by STRING and the GeneMANIA online database. All the DEGs function in metabolic pathways, consistent with the g: Profiler online analysis results. The protein-protein interaction (PPI) networks forecasted by STRING and GeneMANIA were entered into Cytoscape for cytoHubba degree analysis. The hub genes predicted by cytoHubba suggested that fumarate hydratase (FH) might be relevant to DR. qRT-PCR suggested that the expression of FH was higher in DR retinal tissues than in normal control tissues. CONCLUSIONS Multiple bioinformatics analyses verified that FH could be used as a potential diagnostic marker and new therapeutic target of DR.

Background

Diabetes is a serious public health problem worldwide, especially with the increasing obesity rate [1]. Diabetes impairs almost all organ systems, and the progression of diabetic retinopathy (DR) is the main cause of blindness [1,2]. In addition, in late-stage DR, treatment methods such as laser photocoagulation, intraocular steroids, and a PPAR-alpha agonist lipid-lowering drug have not significantly improved quality of life, as previously expected [3,4]. Therefore, early detection and medical intervention are recommended when DR is diagnosed. However, the pathogenesis of DR remains unclear. Current DR markers include the blood levels of glycemia and lipids, which are risk factors, but these biomarkers cannot supply sufficient diagnostic information. Therefore, exploration of new diagnostic biomarkers and novel treatments for DR is needed.

Vascular endothelial growth factor (VEGF) inhibitors were found to be present at a high level in DR patients. Anti-VEGFs provide insight into the treatment of DR [5]. Unfortunately, anti-VEGF treatment does not benefit every DR patient [1]. Therefore, new drug targets would provide further hope for DR patients. With the rapid development of the current gene profiling technique, the identification of differentially expressed mRNAs between patients and healthy controls is easy and convenient. Many mRNA gene profiles have been used to identify different genes associated with many diseases. However, various basic patient parameters such as sex, age, height, weight, and even race have varied among the profile results. Therefore, to date, many differentially expressed genes (DEGs) have been detected by gene profiling, but new biomarkers for DR remain unclear. Recently, the integration of different gene profiles has been successfully used to compensate for these disadvantages.

In our research, by integrating different gene profiles, we removed DEGs that were expressed only in 1 gene profile (“junk information”) [6]. This information might contain some bias due to the different basic parameters such as age, sex, race, and physical condition. Occasionally, the enrichment analysis yielded results that seemed unrelated to the disease. Therefore, we speculated that there was some “junk information” in these DEGs, and attempted to remove the “junk information” from different profiles to identify the truly co-expressed DEGs in DR. In addition, we carried out multiple bioinformatics analyses of the co-expressed DEGs to obtain true information on the core pathogenic genes that could be used as potential diagnostic markers and new therapeutic targets of DR.

Material and Methods

General gene profiling analysis

This study was conducted according to the Declaration of Helsinki and was approved by the Ethics Committee of Jilin University. Raw data on human DR were searched from the Gene Expression Omnibus (GEO) database. The GSE53257 dataset contains 16 tissues (6 DR, 5 normal control and 5 diabetes mellitus without retinopathy tissues), from which we chose 6 DR and 5 normal controls. The GSE60436 dataset contains 9 tissues (3 DR with active neovascularization, 3 DR with no active neovascularization, and 3 normal control tissues), from which we chose 3 DR with active neovascularization and 3 normal controls.

Data processing and integration of 2 profiles

R software (version: R 3.2.3) with the limma package [7] was used to analyze the raw data. The GEO2R method was used to screen the DEGs (version info: Biobase 2.30.0, GEO query 2.40.0, limma 3.26.8). P values less than 0.05 were considered statistically significant. The DEGs with a |fold change| >1.2 were chosen for subsequent analysis. All the up- and downregulated DEGs are shown by volcano plots. All the raw data were displayed in box plots.

Multiple bioinformatics analyses

The Database for Annotation, Visualization and Integrated Discovery (DAVID) online tool (david.ncifcrf.gov, version info: 6.8) was used to predict the biological functions of the DEGs and provide the gene ID for further analysis [8]. Briefly, the DEGs were first entered into the “paste list” section, the “OFFICIAL GENE SYMBOL” was then set as the “select identifier”, and the gene list was selected in the last step. The 2 main biological functions of DAVID are signaling pathway analysis and Gene Ontology (GO) enrichment. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis is a functional analysis that maps gene pathways, and the pathways were verified by g: Profiler online analysis (biit.cs.ut.ee) after gene ID conversion. g: Profiler has been updated with new data from Ensembl and WormBase ParaSite, and the latest version was updated on October 15, 2019. The cellular component (CC), biological process (BP), and molecular function (MF) terms were analyzed by GO analysis in DAVID.

All the DEGs were entered into the STRING online search tool (string-db.org, version info: 11.0) to obtain the physical and functional information of all interacting genes [9,10]. Multiple proteins were selected, and all the DEGs were entered into the “list of names”. The organism selected was Homo sapiens. The basic settings included the meaning of network edges (confidence: line thickness indicates the strength of data support) and active interaction sources (text mining, experiments, databases, co-expression, neighborhood, gene fusion, and co-occurrence). The significant score was selected as greater than 0.4.

The GeneMANIA database () was used to test the results of the GO and pathway enrichment analyses by DAVID. Homo sapiens was set as the main search scope, and all the genes were entered into the “gene list” section. This tool also provides information on the predicted genes relevant to the current DEGs. All the interacting data were entered into Cytoscape software for visualization. Finally, cytoHubba was used to screen the hub genes and signaling pathways. The parameters of cytoHubba were set as follows: Hubba nodes=top 10 nodes ranked by degree, display options=check the first-stage nodes, display the shortest path, and display the expanded subnetwork.

Quantitative reverse transcription-PCR (qRT-PCR) was used to verify the hub genes associated with DR in patient fibrovascular tissues (n=10) and normal retinal tissues (n=10). Ten fibrovascular membranes were surgically removed from 10 eyes of patients with proliferative diabetic retinopathy during pars plana vitrectomy. The mean age was 58.9±4.1 years, there were 4 men and 6 women, and the mean duration of diabetes was 15.7±2.1 years.

Total RNA was reverse-transcribed to cDNA with a PrimeScript RT Reagent Kit (TaKaRa, Japan) according to the manufacturer’s instructions. Primer 5.0 software was used to design primers, and a QuantStudio 7 Flex real-time PCR system (Applied Biosystems, Carlsbad, CA, USA) was used. All samples were normalized to GAPDH. The relative expression levels of each gene were calculated using the 2−ΔΔCt method.

Results

DEGs in DR after raw data processing

In total, 13 085 and 179 DEGs were obtained from the GSE60436 and GSE53257 datasets, respectively. All the normalized data are shown by box plots. By an integrated analysis of these 2 profiles, 37 significant co-expressed DEGs were identified, including 23 downregulated DEGs and 14 upregulated DEGs between DR samples and matched healthy samples. The detailed gene list is shown in Table 1. Volcano plots of DEGs in human DR were generated with microarray technology (Figure 1).

Table 1

The detail gene list of the co-expression genes.

	Gene title	P value	t	B	Log fold change
GSE53257Up	ATP5C1	2.68E-02	2.281599	−5.04323	0.3511132
	PRIM1	1.63E-05	4.773875	1.976374	0.7361976
	BCKDHB	2.05E-02	2.393286	−4.804006	0.379364
	CFHR3	3.06E-02	2.225538	−5.159821	0.3499861
	SLC25A24	3.21E-02	2.205493	−5.200933	0.3802023
	UCHL5	6.40E-03	2.847118	−3.742175	0.4512804
	CLYBL-AS2	2.71E-02	2.277074	−5.052729	0.3590496
	MTHFD1L	2.66E-12	9.193253	17.510878	1.4327238
	AGPAT5	1.39E-04	4.128919	−0.107986	0.6918791
	C10orf2	5.12E-03	2.928792	−3.53677	0.4472831
	KARS	2.37E-02	2.333219	−4.933804	0.3623335
	MTDH	1.82E-05	4.741439	1.868619	0.751102
	SLC25A25	4.19E-02	2.088332	−5.435062	0.3352625
	NOC3L	1.46E-06	5.469513	4.343854	0.8803197
GSE 60436Down	ACYP2	4.81E-04	−3.737387	−1.301455	−0.5895565
	MRPS6	1.42E-06	−5.477544	4.371716	−0.8668793
	PNKD	6.88E-06	−5.025893	2.822394	−0.8232741
	ALAS2	8.31E-06	−4.970738	2.635973	−0.8840675
	ISCU	2.24E-03	−3.223628	−2.76218	−0.498068
	PC	2.93E-14	−10.53169	22.028311	−1.731605
	CYB5R2	5.31E-11	−8.334184	14.514156	−1.2837583
	SLC25A27	1.03E-06	−5.568247	4.687065	−0.8845851
	EHHADH	6.70E-03	−2.829756	−3.785304	−0.4596529
	SLC25A44	1.64E-02	−2.483162	−4.604936	−0.3835724
	SLC25A35	3.03E-02	−2.229506	−5.151647	−0.3458279
	COX17	3.80E-04	−3.813026	−1.075838	−0.6049371
	SPR	5.54E-03	−2.900185	−3.609186	−0.4704401
	ABAT	7.87E-03	2.769408	−3.933723	0.436292
	SLC25A34	5.35E-05	−4.420586	0.81849	−0.6891993
	HMGCS2	1.76E-37	−36.36023	74.834714	−5.9069536
	SLC25A22	1.01E-05	−4.9152	2.448953	−0.7743373
	DUT	4.40E-04	−3.76577	−1.217085	−0.6224354
	ACAT2	1.99E-02	2.405249	−4.777841	0.37544
	ACSL5	7.84E-02	1.796793	−5.970304	0.2769767
	FMC1	1.00E-09	−7.506787	11.577401	−1.1384521
	CPT1B	4.12E-02	−2.095622	−5.420804	−0.32582
	PHYHIPL	1.11E-13	−10.13111	20.697953	−1.8413833
GSE53257Up	ATP5C1	2.37E-04	5.31351	0.8332	0.390735
	PRIM1	2.47E-03	3.88868	−1.4726	0.404072
	BCKDHB	2.68E-03	3.84196	−1.5524	0.410796
	CFHR3	5.89E-03	3.39415	−2.3253	0.377976
	SLC25A24	8.43E-03	3.19308	−2.6752	0.436048
	UCHL5	1.37E-02	2.92311	−3.1445	0.266053
	CLYBL-AS2	1.50E-02	2.87124	−3.2343	0.349242
	MTHFD1L	1.54E-02	2.85703	−3.2588	0.304592
	AGPAT5	1.65E-02	2.81876	−3.3249	0.26755
	C10orf2	2.51E-02	2.58503	−3.7255	0.295845
	KARS	2.57E-02	2.57306	−3.7458	0.387177
	MTDH	3.20E-02	2.44929	−3.9546	0.268641
	SLC25A25	3.47E-02	2.40381	−4.0306	0.47827
	NOC3L	4.24E-02	2.29129	−4.2168	0.288753
GSE53257Down	ACYP2	3.58E-04	−5.04954	0.4276	−0.39137
	MRPS6	6.09E-02	−2.08411	−4.5513	−0.118454
	PNKD	6.31E-04	−4.69708	−0.1303	−0.919907
	ALAS2	3.64E-01	−0.94584	−6.0392	−0.082022
	ISCU	7.40E-04	−4.59982	−0.2873	−0.539286
	PC	8.74E-04	−4.49971	−0.4504	−1.53842
	CYB5R2	1.26E-03	−4.28054	−0.812	−0.32786
	SLC25A27	1.39E-03	−4.22112	−0.911	−1.339729
	EHHADH	2.24E-03	−3.94359	−1.379	−0.758052
	SLC25A44	2.56E-03	−3.86721	−1.5092	−0.587199
	SLC25A35	7.13E-03	−3.28684	−2.512	−1.822979
	COX17	7.20E-03	−3.28116	−2.5219	−0.797026
	SPR	7.98E-03	−3.22403	−2.6213	−0.326796
	ABAT	8.27E-03	−3.20359	−2.6569	−1.740136
	SLC25A34	9.26E-03	−3.14095	−2.766	−1.6249
	HMGCS2	9.46E-03	−3.12871	−2.7872	−0.450437
	SLC25A22	1.28E-02	−2.96085	−3.079	−0.293277
	DUT	1.28E-02	−2.95942	−3.0815	−0.513756
	ACAT2	1.64E-02	−2.82226	−3.3189	−1.127476
	ACSL5	2.25E-02	−2.64639	−3.6209	−0.294467
	FMC1	3.21E-02	−2.4487	−3.9556	−0.266524
	CPT1B	4.21E-02	−2.29536	−4.2102	−0.325261
	PHYHIPL	9.06E-05	−5.95379	1.7721	−0.517559

Figure 1

Comparison of mRNA profile between diabetic retinopathy samples and control samples. (A, B) The volcano plot showed the differentially expressed genes between the diabetic retinopathy and control tissues. The green plot shows decrease and red plot shows increase. (C, D) The box plot showed the distribution of differently expression mRNAs; (E) The volcano plot showed the different co-expression genes after integration of the 2 profiles between the diabetic retinopathy and control tissues.

Bioinformatics analysis results

The co-expressed down- and upregulated genes were predicted by the DAVID online tool. The BP terms of the downregulated DEGs focused on the positive regulation of translation (gene count: 6, P=2.11E-06), transport (gene count: 4, P=3.28E-04), fatty acid beta-oxidation (gene count: 3, P=7.76E-04), metabolic processes (gene count: 3, P=6.25E-03), and cellular iron ion homeostasis (gene count: 3, P=6.25E-03), and the BP terms of the upregulated DEGs focused mainly on transmembrane transport (gene count: 2, P=9.62E-02) (Figure 2). The CC terms of the downregulated focused on the mitochondria (gene count: 11, P=7.26E-08), mitochondrial inner membrane (gene count: 8, P=7.39E-08), and mitochondrial matrix (gene count: 4, P=8.40E-04), and the CC terms of the upregulated DEGs focused on the mitochondrial inner membrane (gene count: 2, P=5.90E-02). The MF terms of the upregulated DEGs focused on heparin binding (gene count: 12, P=7.07E-06), extracellular matrix binding (gene count: 6, P=1.11E-04), integrin binding (gene count: 4, P=3.49E-02), calcium ion binding (gene count: 25, P=3.52E-02), and peptidase activator activity (gene count: 3, P=1.31E-02), and the CC terms of the downregulated DEGs focused on structural constituent of ribosome (gene count: 6, P=1.35E-05) and pyridoxal phosphate binding (gene count: 2, P=8.21E-02).

Figure 2

Gene Ontology bioinformation results.

DAVID showed the BP terms of all co-expressed DEGs focused on translation (gene count: 7, P=1.02E-06), metabolic processes (gene count: 4, P=8.28E-04), transport (gene count: 4, P=1.25E-3), and fatty acid beta-oxidation (gene count: 3, P=1.89E-3), and the CC terms focused on the mitochondrion (gene count: 16, P=8.37E-11), mitochondrial inner membrane (gene count: 9, P=1.36E-07), and mitochondrial matrix (gene count: 5, P=2.09E-04), and the MF terms focused on structural constituent of ribosome (gene count: 7, P=1.07E-05) and catalytic activity (gene count: 4, P=8.26E-04).

The STRING database also showed that the BP terms of all co-expressed DEGs focused on mitochondrial transport (gene count: 9, P=4.03E-07), small molecule metabolic processes (gene count: 18, P=4.03E-07), carboxylic acid metabolic processes (gene count: 11, P=0.13E-03), nucleotide biosynthetic processes (gene count: 0.0033, P=0.33E-2), and monocarboxylic acid metabolic processes (gene count: 7, P=0.004). The CC terms focused on the mitochondrion (gene count: 30, P=7.36E-25), mitochondria (gene count: 22, P=2.17E-17), mitochondrial envelope (gene count: 14, P=1.33E-09), mitochondrial membrane (gene count: 13, P=7.18E-09), and mitochondrial matrix (gene count: 11, P=2.13E-08). The MF terms focused on catalytic activity (gene count: 24, P=0.75E-3), transferase activity (gene count: 6, P=0.75E-3), ATP transmembrane transporter activity (gene count: 2, P=0.59E-2), and ADP transmembrane transporter activity (gene count: 2, P=0.0115).

The enriched pathways revealed by KEGG analysis focused on metabolic pathways (gene count: 16, P=1.02E-08), valine, leucine and isoleucine degradation (gene count: 4, P=7.63E-06), butanoate metabolism (gene count: 4, P=4.99E-05), butanoate metabolism (gene count: 4, P=1.26E-04), and pyruvate metabolism (gene count: 4, P=1.59E-04). g: Profiler revealed valine, leucine, and isoleucine degradation (gene count: 5, P=3.89E-06), butanoate metabolism (gene count: 4, P=2.62E-05), propanoate metabolism (gene count: 4, P=5.88E-05), metabolic pathways (gene count: 13, P=6.83E-05), and pyruvate metabolism (gene count: 3, P=5.88E-03) are the main enrichment pathways. The metabolic pathways accounted for the most genes in both bioinformatics analyses (Figure 3).

Figure 3

Pathways bioinformation results. The gradual color and size mean the P value and gene number.

All the DEGs were entered into the STRING database to obtain the interaction data. The protein-protein interaction (PPI) network exhibited 37 nodes and 24 edges (average node degree: 1.3, PPI enrichment p value: 2.05E-09). The nodes and edges were entered into Cytoscape for further analysis. The degree cytoHubba calculation suggested that the 10 hub genes were fumarate hydratase (FH), EHHADH, ACAT2, ATP5C1, BCKDHB, HMGCS2, PC, SPR, CPT1B, and PRIM1. Among these 10 hub genes, FH accounted for the highest degree score (degree score: 5). The different degree scores of the hub genes are shown as gradual colors in Figure 4.

Figure 4

PPI network of all differentially co-expression genes. (A) PPI of all down- and upregulated co-expression genes. (B) cytoHubba analysis of PPI network with Cytoscape. The gradual color indicates the degree score.

All the co-expressed DEGs, including 23 downregulated DEGs and 14 upregulated DEGs, were entered into the GeneMANIA database for further analysis, including prediction and interaction information analyses. The interaction information on the DEGs analyzed by R software and predicted by GeneMANIA suggested that the mitochondrial matrix acts as the key CC (Figure 5), which was consistent with the DAVID and STRING results. Then, the interaction information was also entered into Cytoscape for cytoHubba analysis. The degree cytoHubba calculation suggested that the 10 hub genes were FH, ACAT2, SLC25A34, PKM, MRPS6, CMC4, MRPS18C, CPT1C, UCHL1, and ACYP2. Among these 10 hub genes, FH had the highest degree score (degree score: 2), which was mostly consistent with the STRING results (Figure 6).

Figure 5

PPI network constructed with GENE Mania. (A) PPI network include all the different co-expression gene and the genes predicted by GENE Mania. (B) PPI network after typed into the co-expression gene. The red node represents the mitochondrial matrix.

Figure 6

PPI network constructed with Cytoscape. (A) The PPI network with the interaction data predicted by GENE Mania. (B) cytoHubba predicted by Cytoscape. The gradual color indicates the degree score.

The primers used in this study are listed in Table 2. qRT-PCR was used to assess the expression of all hub genes analyzed by STRING (DEGs between human DR patient retina tissues and normal control tissues). EHHADH (fold change: 10.23), ACAT2 (fold change: 3.21), HMGCS2 (fold change: 17.02), PC (fold change: 4.31), SPR (fold change: 4.38), and CPT1B (fold change: 3.48) were significantly decreased in human DR patients retina samples, and FH (fold change: 3.68), ATP5C1 (fold change: 3.47), BCKDHB (fold change: 3.74), and PRIM1 (fold change: 18.22) were significantly increased (p<0.05) (Figure 7).

Table 2

Primers used in this study.

Gene name	Forward primer	Reverse primer
FH	GGAGGTGTGACAGAACGCAT	CATCTGCTGCCTTCATTATTGC
EHHADH	AAACTCAGACCCGGTTGAAGA	TTGCAGAGTCTACGGGATTCT
ACAT2	GCGGACCATCATAGGTTCCTT	ACTGGCTTGTCTAACAGGATTCT
ATP5C1	TCACCAGGAGACTAAAGTCCATC	TATTTTGCTGCCGCTACCATT
BCKDHB	GTGCCTTGGATAACTCATTGGC	AGTGCATCTAAAGACTCCACCA
HMGCS2	GACTCCAGTGAAGCGCATTCT	CTGGGAAGTAGACCTCCAGG
PC	ACAGAGGTGAGATTGCCATCC	CACTGCATCTACGTTGTTCTCC
SPR	AAGTGAACAACTACTGGGCAC	AGCATATCACGAGCAGCCTTT
CPT1B	GCGCCCCTTGTTGGATGAT	CCACCATGACTTGAGCACCAG
PRIM1	ATGGAGACGTTTGACCCCAC	CGTAGTTGAGCCAGCGATAGT
GAPDH	CGGACCAATACGACCAAATCCG	AGCCACATCGCTCAGACACC

Figure 7

Validation of qRT-PCR.

Discussion

DR is a serious complication of diabetes that can cause visual impairment in adults. According to previous reports, there are almost 93 million DR patients worldwide, and among these patients, 28 million are at risk of becoming blind [11]. The current treatment strategy is aimed at early diagnosis and treatment, including controlling systemic risk factors, laser therapy, vitrectomy, intravitreal steroids, and anti-VEGF agents [12]. However, the outcomes are not as satisfactory as previously expected. Therefore, novel molecular targets that could be used for early diagnosis and treatment should be researched.

In this study, we used multiple bioinformatics analyses to explore the hub genes in DR. The threshold value in this study was set as 1.2-fold. The threshold value was always set as 2.0-fold between the 2 groups. However, in some situations, there were insufficient co-expressed genes in both profiles to perform enrichment analysis, and the threshold value was set as 1.2-fold [13]. Both the DAVID and STRING online databases suggested that the main BP term of DR is mitochondrial transport, consistent with previous studies [14-16]. Kowluru et al. [14] noted that retinal metabolic abnormalities were a main pathological process during DR that were affected by reactive oxygen species (ROS). ROS can impair mitochondria and mitochondrial DNA, which can lead to abnormal transcription. Then, the electron transport chain is also influenced, and finally, a vicious cycle of free radicals is promoted [14]. Therefore, antioxidants can improve the symptoms of DR patients by inhibiting free radicals. Rodríguez et al. [15] also found that an imbalance between oxidants and antioxidants indeed exists in DR patients. Further research revealed that mitochondrial dysfunction plays an important role during DR [16]. The study noted that oxidative deregulation was decreased in patients with type 2 diabetes mellitus, and membrane fluidity of submitochondrial particles was observed. The multiple bioinformatics analyses described here revealed mitochondrial transport as the main BP term of DR and provided evidence for further research on mitochondrial transport in DR patients.

Profiler and DAVID pathway analyses suggested that metabolic pathways may play a key role during the process of DR. Early treatment by controlling systemic risk factors and regulating the metabolic situation of DR has achieved good outcomes. It had been verified that by controlling the glycemic level, the progression of retinopathy in diabetic patients can be slowed significantly, and the mean risk of retinopathy is reduced by nearly 76% [17]. The UK Prospective Diabetes Study (UKPDS) also noted that metabolic control could have long-term beneficial effects through the “metabolic memory” phenomenon [18]. In addition, metabolic control, such as intensive blood pressure control and angiotensin-converting enzyme inhibitor intervention, has also been proven to have advantages in retinopathy development and progression [19].

After performing GO and signaling pathway analyses, we aimed to verify the hub genes associated with DR. We used CytoHubba in Cytoscape software to analyze the interaction information collected from STRING and GeneMANIA. The 2 online databases revealed that FH can act as a hub gene during DR. The Krebs cycle enzyme FH catalyzes its hydration to malate, and it has been reported to be relevant to tumor suppressors, atypical intradermal smooth muscle neoplasms [20] and nutrient-dependent cytotoxicity [21]. Changing the functional status of FH leads to fumarate accumulation in cells and the activation of HIF-dependent pathways, especially glucose metabolism [22]. Tatsuya et al. [23] first researched the anti-fumarase antibody in diabetic macular edema (DMO) disease, showing that anti-fumarase IgG was related to the visual acuity of patients with DMO disease. In addition, with anti-VEGF treatment, the amount of autoantibody was corrected with visual acuity improvement. Therefore, they inferred, by testing anti-fumarase IgG, that it may become a novel serum biomarker of functional efficiency. Yoshitake et al. [24] also compared serum fumarase IgG in diabetic patients with DMO disease and those without DMO. Fumarase was also relevant to photoreceptor damage and visual dysfunction after administration of the anti-fumarase antibody in animals, and the authors also concluded that autoantibodies against fumarase may participate in the pathological mechanism of DMO disease. Therefore, FH indeed acts as a key gene during the process of DR.

In our study, cholesterol acyltransferase 2 (ACAT2) was also predicted as a hub gene. ACAT is a microsomal protein that can promote cholesterol ester (CE) synthesis. Two subtypes of ACAT exist: ACAT1 and ACAT2. ACAT2 is mainly expressed in the intestine and liver and packages CE into chylomicrons and very low-density lipoproteins [25]. Therefore, ACAT2 can affect hepatic and plasma cholesterol homeostasis and might play roles in the pathogenesis of atherosclerosis [26]. However, the roles of ACAT2 in DR have not been reported. Discovery of the function of ACAT2 in DR might reveal a new pathological mechanism from a new perspective.

Conclusions

Multiple bioinformatics analyses revealed FH as a hub gene during DR. Thus, FH may be used as a diagnostic marker and new drug target. However, some limitations to this study still exist, and the function of FH remains unclear. Further experiments should evaluate the role of FH in metabolic pathways in DR. DAVID and GeneMANIA also revealed ACAT2 as a hub gene. However, this is the first study to suggest that ACAT2 is relevant to DR. Further studies are needed to verify the details of the function of ACAT2 in DR.