Integrated microarray analysis provided novel insights to the pathogenesis of glaucoma.

Glaucoma is characterized as a visual field defect, which is the second most common cause of blindness. The present study performed an integrated analysis of microarray studies of glaucoma derived from Gene Expression Omnibus (GEO). Following the identification of the differentially expressed genes (DEGs) in glaucoma compared with normal control (NC) tissues, the functional annotation, glaucoma‑specific protein‑protein interaction (PPI) network and transcriptional regulatory network constructions were performed. The acute intraocular pressure (IOP) elevation rat models were established and reverse transcription‑quantitative polymerase chain reaction (RT‑qPCR) was performed for DEGs expression confirmation. Three datasets were downloaded from GEO. A total of 97 DEGs, 82 upregulated and 15 downregulated were identified in glaucoma compared with NC groups with false discovery rate <0.05. Response to virus and immune response were two significantly enriched GO terms in glaucoma. Valine, leucine and isoleucine degradation was a significantly enriched pathway of DEGs in glaucoma. According to the PPI network, HDAC1, HBN, UBR4 and PDK1 were hub proteins in glaucoma. FOXD3, HNF‑4 and AP‑1 were the three transcription factors (TFs) derived from top 10 TFs which covered the majority of downstream DEGs in glaucoma. Based on the RT‑qPCR results, the expression levels of 3 DEGs, raftlin, lipid raft linker 1 (RFTN1), PBX homeobox 1 (PBX1), HDAC1 were significantly upregulated and the expression of GEM was significantly downregulated in acute IOP elevation rat model at the first and fifth day. These four DEGs had the same expression pattern with our integrated analysis. Therefore, the current study concluded that 6 DEGs, including HEPH, SELENBP1, RFTN1, ID1, HDAC‑1 and PBX1 and three TFs, including FOXD3, HNF‑4 and AP‑1 may be involved with the pathogenesis of glaucoma. The findings of the current study may improve diagnosis and drug design for glaucoma.

Introduction

As the second most common cause of blindness (1), glaucoma is characterized by a visual field defect induced by progressive loss of retinal ganglion cells (RGCs) and damage to the optic nerve head (ONH) (1,2). Glaucoma may be classified into primary open angle glaucoma (POAG) and closed angle glaucoma according to the angle at the junction between the cornea and iris. POAG is the most common form of glaucoma (3). Glaucoma is a multifactorial disease. Elevated intraocular pressure (IOP) is a major risk factor for glaucoma (4). In addition, old age, family history, racial differences, severe myopia and low diastolic perfusion were also associated with glaucoma (5–8). Various genetic factors, such as vascular endothelial growth factor (VEGF), forkhead box C1 (FOXC1), glutathione S-transferase mu 1 (GSTM1) and cross-linked actin networks (CLANs) has been identified to be associated with glaucoma.

Previous research has shown that lowering IOP is the only clinical treatment demonstrated to slow down the progression of glaucoma (9). Identification of differentially expressed genes (DEGs) in a disease compared with normal control is an approach to find the key genes and pathways associated with the process of diseases and may provide clues for the pathogenesis and identify novel diagnostic and therapeutic strategies. The advance of next generation sequencing technologies allowed for improved identification of DEGs in glaucoma based on previous microarray studies (10–14). Due to the differences in sample size, microarray platforms and analytical techniques among the multiple microarray studies, more accurate profiles of DEGs with larger sample size than an individual microarray can be obtained by integrated analysis of multiple microarrays studies (15).

In order to elucidate the cellular and molecular events associated with glaucoma, the present study identified DEGs associated with glaucoma by integrated microarray analysis of glaucoma. Functional annotation and protein-protein interaction (PPI) network construction were performed to identify glaucoma-associated DEGs and pathways to interpret the biological functions of DEGs. Dysregulation of transcription factors (TFs) has been indicated to be involved with the pathogenesis of various diseases such as glaucoma, the glaucoma-specific transcriptional regulatory network of TFs and DEGs was constructed. In addition, the expression of DEGs was confirmed in the acute IOP elevation rat model. The findings of the current study may provide novel insight in the pathogenesis of glaucoma.

Materials and methods

Eligible gene expression profile of glaucoma

Based on the Gene Expression Omnibus database (GEO; www.ncbi.nlm.nih.gov/geo), the microarray datasets of optic nerve head/optic nerve head astrocytes in glaucoma and normal control (NC) groups were selected. In addition, the selected datasets were the expression profile of whole-genome sequencing and were normalized or original.

Statistical analysis

The present study used the metaMA package in R version 3.2.3 to combine data from multiple microarray datasets and obtained the individual P-values (16). DEGs in glaucoma compared with NC were identified with False discovery rate (FDR) <0.05. The heat-map of DEGs in glaucoma was obtained by using the heatmap.2 function in the Rgplots package (17).

Functional annotation

Functional annotation of DEGs in glaucoma was performed using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway and Gene Ontology (GO) enrichment analyses with the online software GeneCodis3 (genecodis.cnb.csic.es/analysis). FDR <0.05 was considered to be statistically significant.

PPI network constriction

PPI network was used to visualize protein-protein interactions. Using Biological General Repository for Interaction Datasets (BioGRID) (thebiogrid.org/) and Cytoscape version 3.3.0 (www.cytoscape.org/), the PPI network of the top 20 upregulated DEGs and downregulated DEGs in glaucoma were constructed. Nodes were used to represent the proteins and edges were used to represent the interaction between two proteins.

Construction of glaucoma-specific transcriptional regulatory network

The promoter regions of top 20 up and downregulated DEGs in glaucoma were obtained using the University of California Santa Cruz (UCSC) Genome Browser (genome.ucsc.edu). The online tool TRANFAC (18) was used to obtain TFs and motif sequences of binding sites and position weight matrix (PWM). The present study used PWM-scan to obtained the TFs which regulate proteins encoded by DEGs and constructed glaucoma-specific transcriptional regulatory network of top 20 upregulated and downregulated DEGs using Cytoscape version 3.3.0.

Acute elevation of IOP

A total of 9 male Sprague-Dawley rats (age, 8–10 weeks; weight, ~300 g) were housed under normal atmosphere with controlled-temperature (25°C), illumination (12-h light/dark cycle) and with free access to food and water. The 9 rats were divided randomly into control group, model 1 day group and model 5 day group. Except control group, the acute IOP elevation rat model was established by saline perfusion into anterior chamber as previously described (19). IOP was raised by elevating the reservoir of saline 120 cm above the animal's eye for 1 h. IOP was measured using an Icare TonoLab tonometer (ICare Finland Oy, Vantaa, Finland).

RT-qPCR validation of DEGs

After the acute IOP elevation rat model establishment, the rat eyeball tissues in control group (n=3), model 1 day group (n=3) and model 5 day group (n=3) were obtained prior to acute IOP elevation rats model establishment and on the first and fifth day after acute IOP elevation rats model establishment. The total RNA was extracted with TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc., Waltham, MA, USA). cDNA was generated from extracted RNA with SuperScript III Reverse Transcriptase (Invitrogen; Thermo Fisher Scientific, Inc.) incubated at 42°C for 1 h and 72°C for 10 min. In an ABI 7500 real-time PCR system (Applied Biosystems; Thermo Fisher Scientific, Inc.), qPCR was performed with Power SYBR-Green PCR Master mix (Thermo Fisher Scientific, Inc.). Cycling conditions were 95°C for 15 min followed by 40 cycles of 95°C for 10 sec, 55°C for 30 sec and 72°C 32 sec. The relative gene expression of selected DEGs was analyzed using 2−ΔΔCq method (20). The forward and reverse primers (Majorbio Co., Shanghai, China) are presented in Table I. The human 18srRNA was used as endogenous control for mRNA expression in analysis. Quantification of mRNA expression was performed using GraphPad Prism version 5 (GraphPad Inc., La Jolla, CA).

Table I.

Primers used in reverse transcription quantitative polymerase chain reaction.

Gene	Primers (5′-3′)
RFTN1	Forward: agaaggcagagcttcccaat
	Reverse: aggtgttcttcccatccctg
BCKDHB	Forward: atactttaccgggcagcagt
	Reverse: ttttcttgggccatggaagc
CASP1	Forward: acaaagaaggtggcgcattt
	Reverse: aacatcagctccgactctcc
NMI	Forward: atatcgatggggacgcatgt
	Reverse: cctcagacagttcatcggga
PLPP3	Forward: ccttatgtggccgctctcta
	Reverse: tggcccgagaagaaggattt
DCLK1	Forward: ccatcgtgacatcaagccag
	Reverse: aggccatatccagtctctgc
KCTD12	Forward: atggaggagagggagcagta
	Reverse: caaaacgacacccaggacag
HEPH	Forward: accttagctggcaccttgat
	Reverse: tctgtccgtggaagaaagct
ID1	Forward: aacagcaggtgaacgttctg
	Reverse: tcgcgacttcagactcagag
ABCA8	Forward: aaacagcccaccaactccta
	Reverse: ttcatcctgggttgtctgct
MID1	Forward: ttgtgtaaactggttgggcg
	Reverse: cttggcttcttgacgggatg
PBX1	Forward: agccgaattgcaggtctttc
	Reverse: tctcctcttccagccctttg
SEPP1	Forward: ggtttgccctactccttcct
	Reverse: ttgtcatggtgcttgtggtg
FUT8	Forward: gaaaattcacttcggggcgt
	Reverse: tctggttgtgggcattttgg
NBN	Forward: cagagctggcagttcaagtg
	Reverse: tcctcactgctgtcctgaag
SELENBP1	Forward: atacatgtgtgggactggca
	Reverse: gtagaagcgctggatgttgg
HDAC1	Forward: tacgacggggatgttggaaa
	Reverse: ttggtcatctcctcagcgtt
AMIGO2	Forward: tgccatgttccaggagctaa
	Reverse: agatcagccagcttgaacct
GEM	Forward: tctgcggtggaagagttgat
	Reverse: tctttgggtcaatctgggct
CYFIP2	Forward: acttcctccccaactactgc
	Reverse: gcggtaggagctgtagatgt

RFTN1, raftlin, lipid raft linker 1; PBX1, PBX homeobox 1; ID1, inhibitor of differentiation 1; SELENBP1, selenium binding protein 1.

Results

Identification of DEGs in glaucoma

Three datasets, including GSE45570, GSE9944 and GSE2378 were used in the current study (Table II). A total of 97 DEGs (82 upregulated and 15 downregulated were identified in glaucoma samples compared with the NC group with FDR <0.05. The heatmap of DEGs in glaucoma was illustrated in Fig. 1 and the top 20 DEGs in glaucoma were displayed in Table III.

Table II.

mRNA expression datasets used in this study.

GEO accession no.	Platform	Normal:Case	Year	Tissue
GSE45570	GPL5175	6:6	2013	Optic nerve head
GSE9944	GPL8300	6:13	2008	Optic nerve head astrocytes
GSE2378	GPL8300	3:7	2005	Optic nerve head astrocytes

Figure 1.

Heatmap of differentially expressed genes in glaucoma samples compared with normal control.

Table III.

Top 20 differentially expressed genes in glaucoma compared with normal control.

ID	Gene name	Combined ES	P-value	FDR	Regulation
23180	RFTN1	1.97057	1.18×10−6	0.001965	Up
594	BCKDHB	1.863234	1.38×10−6	0.001965	Up
834	CASP1	1.901172	4.48×10−7	0.001965	Up
9111	NMI	1.813946	1.10×10−6	0.001965	Up
347902	AMIGO2	−1.86914	2.38×10−6	0.002707	Down
8613	PLPP3	1.69169	3.85×10−6	0.002979	Up
9201	DCLK1	1.84913	4.72×10−6	0.002979	Up
115207	KCTD12	1.796888	5.38×10−6	0.003057	Up
2669	GEM	−1.66006	6.53×10−6	0.003149	Down
9843	HEPH	1.663202	7.21×10−6	0.003149	Up
3397	ID1	1.677848	8.45×10−6	0.003428	Up
10351	ABCA8	1.753583	1.68×10−5	0.006363	Up
4281	MID1	1.770553	1.81×10−5	0.006411	Up
5087	PBX1	1.616593	2.08×10−5	0.006948	Up
6414	SEPP1	1.711474	2.27×10−5	0.007171	Up
2530	FUT8	1.5493	2.47×10−5	0.007393	Up
3065	HDAC1	1.642544	3.09×10−5	0.008781	Up
26999	CYFIP2	−1.5543	4.16×10−5	0.010273	Down
4683	NBN	1.51878	3.96×10−5	0.010273	Up
8991	SELENBP1	1.577832	3.99×10−5	0.010273	Up

RFTN1, raftlin, lipid raft linker 1; FDR, false discovery rate; ES, effect size; PBX1, PBX homeobox 1; SELENBP1, selenium binding protein 1.

According to the GO enrichment analysis, response to virus (FDR=7.11×10−5), immune response (FDR=0.00252621), protein binding (FDR=2.59×10−10), transcription factor binding (FDR=0.000806942), nucleus (FDR=1.24×10−6), extracellular space (FDR=0.000618804) were significantly enriched GO terms of DEGs in glaucoma. The top 20 significantly enriched GO terms including biological process, cellular component and molecular function of DEGs in glaucoma were displayed in Fig. 2. After the KEGG enrichment analysis, Valine, leucine and isoleucine degradation (FDR=0.0213357), Cytokine-cytokine receptor interaction (FDR=0.0220177), Glycolysis/gluconeogenesis (FDR=0.0293935) were most significantly enriched pathways of DEGs in glaucoma. The significantly enriched KEGG pathways of DEGs in glaucoma were displayed in Fig. 3 and Table IV.

Figure 2.

Top 20 significantly enriched Gene Ontology terms of differentially expressed genes in glaucoma compared to normal control.

Figure 3.

Top 20 significantly enriched Kyoto Encyclopedia of Genes and Genomes pathways of differentially expressed genes in glaucoma compared to normal control. BP, biological processes; MF, molecular function, CC, cellular components; FDR, false discovery rate.

Table IV.

Significantly enriched Kyoto Encyclopedia of Genes and Genomes pathways of differentially expressed genes in glaucoma compared with normal control.

ID	Function	Count	P-value	FDR	Genes
00280	Valine, leucine and isoleucine degradation	3	0.0002452	0.0213357	ALDH7A1, BCKDHB, ALDH2
04060	Cytokine-cytokine receptor interaction	5	0.0007592	0.0220177	CXCL6, IL15, IFNAR2, CXCL12, TNFSF10
00010	Glycolysis/gluconeogenesis	3	0.0006757	0.0293935	ALDH7A1, PFKP, ALDH2
00410	β-alanine metabolism	2	0.0023259	0.0337254	ALDH7A1, ALDH2
00053	Ascorbate and aldarate metabolism	2	0.0023259	0.0337254	ALDH7A1, ALDH2
05323	Rheumatoid arthritis	3	0.0016298	0.0354476	CXCL6, IL15, CXCL12
00340	Histidine metabolism	2	0.0028901	0.0359192	ALDH7A1, ALDH2
00640	Propanoate metabolism	2	0.0035119	0.0381914	ALDH7A1, ALDH2
04672	Intestinal immune network for IgA production	2	0.0065563	0.043877	IL15, CXCL12
00300	Lysine biosynthesis	1	0.0082213	0.0447033	ALDH7A1
04120	Ubiquitin mediated proteolysis	3	0.0062022	0.0449656	CUL4B, UBE2L6, MID1
00071	Fatty acid metabolism	2	0.0057135	0.0451886	ALDH7A1, ALDH2
00380	Tryptophan metabolism	2	0.0057135	0.0451886	ALDH7A1, ALDH2
00561	Glycerolipid metabolism	2	0.0074521	0.0463094	ALDH7A1, ALDH2
00310	Lysine degradation	2	0.0080782	0.0468536	ALDH7A1, ALDH2
00330	Arginine and proline metabolism	2	0.0097429	0.0470908	ALDH7A1, ALDH2
05130	Pathogenic Escherichia coli infection	2	0.0097429	0.0470908	NCK1, LY96

FDR, false discovery rate.

Glaucoma-specific PPI network

A PPI network (Fig. 4) of the top 20 upregulated DEGs in glaucoma was constructed, which consisted of 813 nodes and 865 edges. The hub proteins were HDAC1 (degree=456), HBN (degree=75) and MID1 (degree=61). The PPI network (Fig. 5) of the downregulated DEGs in glaucoma was constructed, which consisted of 372 nodes and 393 edges. The hub proteins were UBR4 (degree=56), PDK1 (degree=45) and PHKG2 (degree=41).

Figure 4.

Glaucoma-specific protein-protein interaction network of top 20 upregulated DEGs. The red ellipses indicate the proteins encoded by upregulated DEGs and the purple ellipses indicate other proteins. DEGs, differentially expressed genes.

Figure 5.

Glaucoma-specific protein-protein interaction network of top 20 downregulated DEGs. The green ellipses indicate the proteins encoded by downregulated DEGs and the purple ellipses indicate the other proteins. DEGs, differentially expressed genes.

Glaucoma-specific transcriptional regulatory network

The glaucoma-specific transcriptional regulatory network of top 20 upregulated DEGs (Fig. 6) was consisted of 57 nodes and 141 edges. The glaucoma-specific transcriptional regulatory network of downregulated DEGs (Fig. 7) consisted of 53 nodes and 103 edges. The top 10 TFs covering the majority of downstream DEGs in glaucoma were presented in Table V.

Figure 6.

Transcriptional regulatory network of top 20 upregulated DEGs in glaucoma. The red ellipses indicate the upregulated DEGs and the blue rhombus indicate the TFs. The lines represent the TF-DEG pairs. DEGs, differentially expressed genes; TF, transcription factor.

Figure 7.

Transcriptional regulatory network of downregulated DEGs in glaucoma. The green ellipses represented the downregulated DEGs and the blue rhombus indicated the TFs. The lines represent the TF-DEG pairs. DEGs, differentially expressed genes; TF, transcription factor.

Table V.

Top 10 TFs covering the majority of downstream DEGs in glaucoma.

TF	Count	DEGs
Pax-4	50	SERPING1, SV2A, EFHD1, NMI, DCLK1, CRELD1, CRELD1, AMIGO2, NBN, SERPING1, NBN, NBN, SLC12A2, NBN, FUT8, SELENBP1, SERPING1, SLC25A1, EFHD1, FUT8, LY96, EFHD1, HEPH, ID1, SERPING1, CYFIP2, UBR4, ABCA8, HDAC1, GEM, ID1, SELENBP1, PHKG2, SERPING1, EFHD1, CRELD1, CRELD1, ID1, LY96, SERPING1, ID1, LY96, HDAC1, EFHD1, SELENBP1, LY96, DMD, MID1, HDAC1
1-Oct	44	DMD, DMD, MID1, DMD, NBN, DMD, MID1, PDK1, BCKDHB, BCKDHB, NBN, DMD, BCKDHB, DMD, MID1, DMD, KCTD12, PLPP3, NMI, C1R, DMD, PBX1, SERPING1, DMD, UBR4, SERPING1, MID1, DCLK1, KCTD12, ABCA8, PDK1, SERPING1, DCLK1, PBX1, KCTD12, DMD, PHKG2, PBX1, NBN, BCKDHB, MID1, CASP1, NBN, DMD
Nkx2-5	38	HEPH, NBN, SV2A, NBN, NBN, SERPING1, NBN, SERPING1, ABCA8, EFHD1, DMD, HEPH, SELENBP1, LY96, DCLK1, SELENBP1, PHKG2, SV2A, HDAC1, SLC12A2, BCKDHB, EFHD1, DMD, BNIP3, EFHD1, SERPING1, AMIGO2, RFTN1, SERPING1, ABCA8, ABCA8, SELENBP1, UBR4, EFHD1, LY96, KCTD12, ABCA8, KCTD12
COMP1	26	SLC12A2, MID1, MID1, SLC25A1, MID1, UBR4, DMD, PBX1, HEPH, FUT8, DMD, MID1, NBN, DMD, NBN, MID1, NBN, DCLK1, NBN, PHKG2, DMD, MID1, PBX1, SELENBP1, SLC25A1, DCLK1
FOXD3	23	DMD, DMD, SLC25A1, PBX1, PLPP3, SERPING1, PLPP3, DMD, DCLK1, DCLK1, SERPING1, PLPP3, DMD, PDK1, DMD, DMD, SERPING1, LY96, SLC25A1, CYFIP2, UBR4, DMD, PDK1
HNF-1	23	ABCA8, DMD, DMD, DMD, RFTN1, PHKG2, HEPH, DMD, SV2A, DMD, CYFIP2, DMD, AMIGO2, MID1, SELENBP1, FUT8, CASP1, NBN, NBN, KCTD12, DCLK1, DCLK1, DMD
HNF-4	20	MID1, GEM, PDK1, BNIP3, BCKDHB, CYFIP2, SLC25A1, UBR4, MID1, SLC25A1, BNIP3, DMD, MID1, SLC25A1, KCTD12, CYFIP2, SLC12A2, PDK1, SLC25A1, CASP1
Evi-1	18	HDAC1, DMD, MID1, BCKDHB, FUT8, EFHD1, LY96, DMD, DMD, BCKDHB, DMD, MID1, PBX1, PBX1, SERPING1, SELENBP1, PHKG2, FUT8
AP-1	15	NBN, HEPH, HEPH, PHKG2, MID1, HDAC1, RFTN1, PHKG2, BCKDHB, DCLK1, CYFIP2, DCLK1, FUT8, PHKG2, BNIP3
HNF-3β	12	PBX1, SERPING1, LY96, NMI, SELENBP1, DMD, SERPING1, DCLK1, ABCA8, DMD, CASP1, DMD

TFs, transcription factor; DEGs, differentially expressed genes; PBX1, PBX homeobox 1; RFTN1, raftlin, lipid raft linker 1; ID1, inhibitor of differentiation 1.

Acute elevation of IOP and RT-qPCR validation of DEGs

IOP was significantly increased in model 1 day group and model 5 day group after the acute IOP elevation rats model establishment. In order to verify the expression of integrated analysis, the top 20 DEGs in glaucoma were investigated. Based on the RT-qPCR (Fig. 8), the expression of 3 DEGs [raftlin, lipid raft linkerm 1 (RFTN1), PBX homeobox 1 (PBX1), HDAC1] were significantly upregulated and the expression of GEM was significantly downregulated in the acute IOP elevation rat model on the first and fifth day. These four DEGs presented the same pattern in the integrated analysis in the present study.

Figure 8.

Reverse transcription-quantitative polymerase chain reaction results of the top 20 differentially expressed genes in acute intraocular pressure elevation rat models. *P<0.05 and **P<0.01 vs. control group.

Discussion

In order to reveal the cellular and molecular events associated the pathogenesis of glaucoma, the present study identified several genes and TFs that were involved with glaucoma based on the DEGs in glaucoma compared with the NC group by using an integrated microarray analysis.

The acute IOP elevation rat model was established to validate the expression of DEGs in glaucoma using RT-qPCR. The acute IOP elevation model has been previously used to simulate acute angle-closure glaucoma (21,22); however, the association between acute IOP elevation and chronic glaucoma remains to be fully elucidated. Acute ischemic insult and non-ischemic injuries may be induced by the acute IOP elevation model and lead to significant loss of RGCs after 3 days of acute IOP elevation (23). Alterations in the morphology of astrocytes and influence on axons and cellular responses within the ONH may also be induced by acute IOP elevation (24,25). Additionally, Morrison et al previously reported that gene expression profiles within the ONH induced by acute IOP elevation were similar to those observed following chronic IOP elevation (26,27). Due to the similarity between the acute IOP elevation model and glaucoma, the gene expression profiles in acute IOP model may identify the key genes in glaucoma.

A previous study has indicated that the ONH is repeatedly exposed to cycles of relatively minor injuries in the early stage of glaucoma (28) and these minor injuries may be reversed by treatment with healthy RGCs (28). However, repeated insults or an intrinsic impairment may overwhelm the capacity of RGCs to recover, and the cell death program may be induced (29). Although 1 day of acute IOP elevation is insufficient to induce evident damage on the ONH, the present study determined that the process of RGC recovery is activated after 1 day of acute IOP elevation which may be associated with early-stage glaucoma. Therefore, the expression of selected DEGs after 1 day and 5 days was validated.

As the major risk factor of glaucoma, elevated IOP may induce the optic nerve axonal compression at the lamina cribrosa, block the axoplasmic flow, cause interference in retrograde neurotrophin transport to RGCs, and drive apoptosis. This may lead to disruption of the balance between aqueous humor production and outflow from the anterior chamber (30). Elements regulating IOP and aqueous humor may be involved with the pathogenesis of glaucoma.

As a transporter expressed by ciliary body epithelia, RFTN1 was believed to contribute to the pathogenesis of POAG by changing the composition of aqueous humor (31). In addition, previous studies have reported that raftlin, RFTN1 rs690037 was associated with three glaucoma-associated optic disc parameters including optic cup area, vertical cup-to-disc ratio and central corneal thickness (32,33). Combination of RFTN1 and other optic disc parameter-associated genes, such as ATOH7 single nucleotide polymorphisms conferred risk for POAG (32). In the integrated analysis of the present study, RFTN1 was the top upregulated DEG in glaucoma compared with NC and it was also upregulated in the acute IOP elevation rat model, which provided evidence for the potential role of RFTN1 in glaucoma. In order to determine the specific role of RFTN1 in glaucoma, further research is required.

Hephaestin (HEPH) and selenium binding protein 1 (SELENBP1) were previously reported to be associated with elevated intraocular pressure (34) and upregulated HEPH has been previously detected in glaucoma (35). In the current study, HEPH and SELENBP1 were also upregulated DEGs in glaucoma and the acute IOP elevation rat model. Therefore, it is possible that HEPH and SELENBP1 may be involved with glaucoma by regulating IOP.

However, lowering IOP is not always an effective treatment for patients with glaucoma, which suggested that other mechanisms of glaucoma require further research.

Glaucoma is the most common optic neuropathy and RGCs apoptosis is a major hallmark of glaucoma (36). The present study identified several DEGs in glaucoma which were associated with RGCs damage and apoptosis.

Retinal injury, including light damage and N-Methyl-D-aspartic acid (NMDA) treatment have been reported to induce an increase in bone morphogenetic protein (BMP) expression and BMP-Smad1/5/8 signaling in inner retinal cells, which suggested that BMP-Smad1/5/8 signaling had a neuroprotective effect on RGCs after damage (37). As a target of BMP-Smad1/5/8, inhibitor of differentiation 1 (ID1) was upregulated in response to retinal injury (37). According to the present study, ID1 was also upregulated in the acute IOP elevation rat model and glaucoma which is a RGCs disease. The current study concluded that ID1 was involved with RGC damage, which may be a potential target for neuroprotective therapies in glaucoma and other RGC diseases.

Nuclear atrophy is one of the early events of RGCs death and previous studies have reported that Class I histone deacetylases (HDACs) 1, 2 and 3 have key roles in nuclear atrophy and apoptosis of RGCs (36). In addition, inhibition of HDACs activity in the RGCs nucleus may prevent RGCs death. Treatment of purified rat RGCs with inhibitors including valproic acid (VPA) and sodium butyrate (SB), prevented histone deacetylation and protected RGCs from death in vitro (38). MS-275, an inhibitor of HDAC1/HDAC3 was identified to have a protective effect on the loss of RGCs and promoted RGC survival following optic nerve injury (39). According to the PPI network of the top 20 upregulated DEGs in glaucoma, HDAC1 is a hub protein which was upregulated in the acute IOP elevation rat model as well. The study concluded that inhibitors of HDACs, particularly the HDAC1 inhibitor may have a protective role of RGCs in glaucoma (39). Further investigation is required to identify the specific roles of each individual HDAC in injured RGCs and determine which HDAC may be the appropriate target for therapeutic inhibition in optic neuropathies such as glaucoma.

According to the GO and KEGG enrichment analysis, various changes occurred in glaucoma compared with NC and the majority of enriched pathways in glaucoma were metabolic-associated pathways. Glycolysis/gluconeogenesis and valine, leucine and isoleucine degradation were two significantly enriched pathways in glaucoma which were consistent with previous studies (40–42) and may have important roles in glaucoma. These metabolic-associated pathways may involve with the process of glaucoma and these abnormal metabolites may make a contribution in early diagnosis and treatment of glaucoma.

The present study used the glaucoma-specific transcriptional regulatory network, and identified three TFs which may be involved with the process of glaucoma. Ocular malformation is another risk factor of glaucoma, which may be induced by dysplasia of the anterior eye (43). The structures of anterior eye are closely associated with the migratory neural crest cells. The FOXD3 encodes a forkhead TF, which has a crucial role in neural crest specification in vertebrates; therefore, FOXD3 may be involved in multiple types of eye diseases, such as glaucoma (44). Previous studies have indicated that FOXC1 has a role in various glaucoma phenotypes (45,46). Additionally, another upregulated DEG in glaucoma and the acute IOP elevation rat model, PBX1 may impair the activity of FOXC1 in a filamin A-mediated manner. Considering the present findings and previous studies, it is possible that PBX1 may be involved with the process of glaucoma by interacting with FOXC1. HNF-4 and AP-1 were two TFs identified from top 10 TFs covering the majority of downstream DEGs which have a role in cellular regulation of the reactive ONH astrocytes (11). Therefore, the present study concluded that HNF-4 and AP-1 may be closely associated with glaucoma.

Additionally, to the best of our knowledge this is the first study to report the association between glaucoma and other 6 DEGs, including NMI, KCTD12, ABCA8, MID1, FUTB, and GEM, which were also confirmed by RT-qPCR in the acute IOP elevation rat model; however, their potential roles in glaucoma require further investigation.

In conclusion, HEPH, SELENBP1 and RFTN1 may be involved with the pathogenesis of glaucoma by regulating IOP or glaucoma-associated optic disc parameters. ID1 and HADC1 may be associated with glaucoma by the influence on RGCs. Interaction between PBX1 and FOXC1 may have a role in glaucoma. Three TFs, including FOXD3, HNF-4 and AP-1 may also act as regulators of glaucoma. The findings of the present study may contribute to developing novel and more effective approaches of diagnosis and drug design for glaucoma. Future investigations of glaucoma-associated genes, require the chronic IOP elevation model and the specific role of these glaucoma-associated DEGs and TFs in glaucoma necessary in further research.