This article is mentioned in:

Introduction

Glaucoma is a widely known, multi-factorial disease, which may result in apoptosis of retinal ganglion cells. According to the World Health Organization, glaucoma is the second principal cause of blindness and the most common cause of irreversible blindness in the world (1,2). In glaucoma, the anterior and posterior segments of the eye are affected, and serious damage may be detected in the trabecular meshwork (3). Oxidative stress is considered to be responsible for the molecular damage in the anterior chamber. Primary open-angle glaucoma (POAG) is the most common type of glaucoma, accounting for 60–70% all glaucoma (4). A candidate protein that may be associated with POAG is myocilin (MYOC), encoded by the MYOC gene. MYOC mutations are common in patients with POAG with high levels of intraocular pressure (IOP) (5,6). Additionally, mutations in optineurin were identified in patients with POAG (7). Previous studies suggested that an abnormal expression of serine/threonine-protein kinase TBK1 is a cause of normal-tension glaucoma (8–10). Furthermore, a previous study suggested that the calcium load-activated calcium channel was involved in glaucoma and that cyclin-dependent kinase 4 inhibitor B antisense RNA 1 was upregulated in the retina of a rat model of glaucoma (11).

However, even substantial decreases in IOP are not able to prevent the development and progression of glaucoma in a number of clinical cases (12). Glaucoma-associated cell death is primarily caused by apoptosis, which is triggered by oxidative stress via mitochondrial damage, inflammation, endothelial dysregulation and dysfunction, and hypoxia (13). In general, glaucoma is not preventable; however, the vast majority of patients may maintain useful visual function for life if they have early detection and appropriate treatment (14). Therefore, for the prevention of glaucoma, emphasis must be placed on early detection, and early diagnosis and treatment.

The rapid development and application of high-throughput sequencing technology has provided a comprehensive and rapid analytical method for the study of the pathogenesis of glaucoma, and provide novel ideas for the future treatment of glaucoma (15). The present study aimed to analyze high-throughput transcriptome data from tissue samples of patients with glaucoma and a normal control group. The data was used in bioinformatics analyses to identify key transcription factors (TFs) associated with glaucoma, to examine the pathogenesis of glaucoma and provide a basis for the diagnosis of glaucoma and drug development.

Materials and methods

Microarray expression profiling in Gene Expression Omnibus (GEO)

The GEO is the largest database of high-throughput gene expression data that was developed and is maintained by the National Center for Biotechnology Information (16). The GEO was searched to obtain gene expression profiling studies of glaucoma subjects. The following key search terms were used: [‘glaucoma’ (Medical subject headings Terms) OR ‘glaucoma’ (All Fields)] AND ‘Homo sapiens’ (porgn) AND ‘gse’ (Filter). The selection criteria were as follows: i) The selected dataset must include genome-wide mRNA transcriptome data; ii) the data was obtained from the trabecular meshwork tissue samples of glaucoma and normal control trabecular meshwork tissue samples; and iii) normalized and raw datasets were considered. Following selection, two sets of GSE27276 (17) and GSE4316 (18) glaucoma mRNA data were obtained (19,20).

Identification of differentially expressed genes (DEGs) in glaucoma compared with normal controls

Background correction was performed on the raw data. The normalization was performed using the Linear Models for Microarray (Limma version 3.30.13) Data package in R (21). Subsequently, two-tailed Student's t-tests were performed to calculate individual P-values. Stouffer's test was used to merge individual P-values, and multiple comparison correction was performed using the Benjamini and Hochberg method to obtain the false discovery rate (FDR) (22). Genes with FDR<0.001 were selected as DEGs. Finally, the DEGs in glaucoma vs. normal were identified.

Functional annotation of DEGs

Gene Ontology (GO) (23) and Kyoto Encyclopedia of Genes and Genomes (KEGG) (24) pathway enrichment analysis were performed to detect the biological functions and potential pathways associated with DEGs using GeneCoDis3 (http://genecodis.cnb.csic.es/analysis) as previously described (19). The GO functions of the DEGs were determined according to the three categories of: ‘Biological process’; ‘molecular functions’; and ‘cellular component’. Pathway enrichment analysis was based on the KEGG database, as previously described (25).

Protein-protein interaction (PPI) network construction

In order to identify candidate genes involved in the formation of glaucoma, PPI networks of significant DEGs were constructed, according to the data from Biological General Repository for Interaction Datasets (BioGRID; http://thebiogrid.org/). Based on the existing data of protein interaction in the BioGRID database, Cytoscape (www.cytoscape.org/; (version 3.5.0) was used to search the top 100 upregulated and downregulated mRNAs. The PPI network interaction map was generated subsequent to the exclusion of genes that were not differentially expressed.

Screening for TFs of the top 20 DEGs and construction of TF regulatory network

For the top 20 DEGs, the 2 kb upstream promoter regions were downloaded from the University of California Santa Cruz (UCSC) Genome Browser website (genome.ucsc.edu). The TRANScription FACtor (TRANSFAC) website match tool (gene-regulation.com/pub/databases.html) was subsequently used to analyze TFs capable of binding to the promoter region of the DEGs. TFs that exhibited altered expression in glaucoma with FDR<0.001 were selected. Position Weight Matrix scanning was used to scan the human genome sequence to obtain the protein-coding genes that were regulated by the differentially expressed TFs. Following removal of redundant information, the glaucoma-specific transcriptional regulatory network was constructed using Cytoscape software.

In silico validation of DEGs using GEO

The GEO database (GSE9944) was used to validate the expression of selected glaucoma DEGs. The expression levels of these genes were compared between the glaucoma cases and the normal group. The expression of five genes, cAMP responsive element binding protein 1 (CREB1), fibronectin 1 (FN1), keratin 19 (KRT19), lipocalin 2 (LCN2) and paired box 6 (PAX6) was detected, with the difference of expression levels presented as box-plots.

Results

Differential expression analysis of genes in glaucoma

The probes corresponding to multiple genes were removed, and the average gene expression to which multiple probes corresponded with was calculated. Finally, the intersection of 15,757 genes was obtained.

A total of two gene expression microarray datasets (GSE27276 and GSE4316) were used for the analysis. Compared with the normal controls, 1,935 DEGs in glaucoma were obtained (P<0.05); among these, 951 genes were upregulated and 984 genes were downregulated. The top 40 most significantly up- or downregulated genes are summarized in Table I. Among which, PAX6 (26), LCN2 (27), and MAOA (28) were downregulated and were associated with glaucoma. The DEGs were screened for clustering analysis. The heatmap produced by cluster analysis of the two sets of cDNA microarray data is presented in Fig. 1.

Figure 1. Heatmap of genes that were significantly up- or downregulated in glaucoma compared to normal controls. Row and column represent DEGs and tissue samples, respectively. The color scale indicated the expression of DEGs. The expression data was log transformed and standardized to zero mean and unit variance. Red and green color indicate up- and down-regulation, respectively. DEGS, differentially expressed genes.

Table I. Top 40 differentially expressed mRNAs.

Table I.

Top 40 differentially expressed mRNAs.

A, Upregulated
ID	Symbol	Combined.ES	P-value	FDR
3043	HBB	5.16018337	0	0
9415	FADS2	2.35332790	1.40×10−9	3.68×10−6
56172	ANKH	2.33974344	7.72×10−9	1.11×10−5
3113	HLA-DPA1	2.22895504	1.17×10−8	1.41×10−5
5538	PPT1	2.27584322	1.72×10−8	1.69×10−5
83643	CCDC3	2.09336726	2.56×10−8	2.32×10−5
81618	ITM2C	2.10210833	5.11×10−8	3.81×10−5
8857	FCGBP	2.00251174	1.07×10−7	6.25×10−5
4166	CHST6	1.94047822	1.48×10−7	7.07×10−5
4753	NELL2	1.97462732	2.45×10−7	1.07×10−4
8718	TNFRSF25	1.80276729	8.11×10−7	2.84×10−4
81552	VOPP1	1.82118066	8.70×10−7	2.91×10−4
2191	FAP	1.74786612	1.01×10−6	3.09×10−4
187	APLNR	1.74800187	1.02×10−6	3.09×10−4
51196	PLCE1	1.77349072	1.23×10−6	3.61×10−4
23676	SMPX	1.76299253	1.37×10−6	3.78×10−4
51226	COPZ2	1.72176939	1.72×10−6	4.24×10−4
1290	COL5A2	1.71261435	1.78×10−6	4.31×10−4
57559	STAMBPL1	1.69907540	1.81×10−6	4.31×10−4
2331	FMOD	1.75952553	1.83×10−6	4.31×10−4
B, Downregulated
ID	Symbol	Combined.ES	P-value	FDR
3880	KRT19	−2.54377891	3.23×10−10	1.70×10−6
116039	OSR2	−2.44873456	7.08×10−10	2.23×10−6
5648	MASP1	−2.39197514	2.01×10−9	4.52×10−6
79845	RNF122	−2.31970767	3.00×10−9	5.91×10−6
3934	LCN2	−2.20821475	8.44×10−9	1.11×10−6
57801	HES4	−2.22824712	1.65×10−8	1.69×10−5
10232	MSLN	−2.09546485	2.65×10−8	2.32×10−5
9245	GCNT3	−2.11365461	4.07×10−8	3.38×10−5
84525	HOPX	−2.01200085	5.57×10−8	3.81×10−5
4128	MAOA	−2.03058309	7.95×10−8	5.22×10−5
5080	PAX6	−2.09500899	8.92×10−8	5.62×10−5
64073	C19orf33	−2.02417662	1.26×10−7	6.86×10−5
3866	KRT15	−1.96990194	1.40×10−7	6.92×10−5
7148	TNXB	−1.92147744	1.40×10−7	6.92×10−5
163732	CITED4	−1.94673475	1.94×10−7	9.01×10−5
1675	CFD	−1.97823285	2.68×10−7	1.14×10−4
6510	SLC1A5	−1.95192053	2.77×10−7	1.15×10−4
1638	DCT	−1.84682776	3.22×10−7	1.30×10−4
3861	KRT14	−1.77674842	7.12×10−7	2.61×10−4
10053	AP1M2	−1.77498212	8.52×10−7	2.92×10−4

[i] FDR, false discovery rate.

Functional annotation

In Fig. 2, GO enrichment demonstrated that the DEGs were significantly enriched the ‘Biological processes’ categories: ‘eye development’ (FDR=4.15×10−3); ‘visual perception’ (FDR=7.13×10−3); ‘negative regulation of insulin receptor signaling pathway’ (FDR=1.47×10−2); the ‘Cellular components’ categories: ‘membrane’ (FDR=3.41×10−52); ‘endoplasmic reticulum’ (FDR=3.76×10−12); ‘cytoplasm’ (FDR=1.73×10−63); and ‘Molecular functions’ categories: ‘nucleotide binding’ (FDR=8.50×10−16); ‘hydrolase activity’ (FDR=8.50×10−16); and ‘protein binding’ (FDR=6.15×10−60).

Figure 2. Gene ontology functional enrichment of differentially expressed genes (FDR<0.05). (A) BP functions. (B) CC functions. (C) MF functions. FDR, false discovery rate; BP, Biological process; CC, Cellular component; MF, Molecular functions.

Furthermore, as presented in Fig. 3, the results of the KEGG pathway enrichment demonstrate that DEGs were enriched in ‘Staphylococcus aureus infection’ (FDR=6.15×10−60); ‘Pathways in cancer’ (FDR=8.66×10−5); ‘Systemic lupus erythematosus’; and ‘Type I diabetes mellitus’ (FDR=0.84×10−3).

Figure 3. Kyoto Encyclopedia of Genes and Genomes analysis of differentially expressed genes in glaucoma. FDR, false discovery rate.

PPI network construction

Following removal of the non-DEGs, the PPI network was established. The results are presented in Fig. 4. The network consisted of nodes and edges. The nodes in the network represent the proteins and the edges represent the interactions between them. There were 290 nodes and 290 edges identified. Among them, the genes with higher degrees were: FN1 (degree=92); KRT15 (degree=19); major histocompatibility complex, class II, DP α 1 (HLA-DPA1; degree=12); CREB1 (degree=11); KRT19 (degree=9); histone cluster 2 H2B family member e (degree=9); cell division cycle associated 7 like (degree=7); PAX6 (degree=7); FOS like 2, AP-1 transcription factor subunit (degree=7); minichromosome maintenance complex component 2 (degree=6); palmitoyl-protein thioesterase 1 (degree=6); and lysine demethylase 5B (degree=6). The hub proteins were FN1 (degree=92), KRT15 (degree=19) and HLA-DPA1 (degree=12).

Figure 4. Protein-protein interaction network of top 100 significantly DEGs in glaucoma. The blue ellipses represent the proteins encoded by downregulated DEGs and the red ellipses represent the proteins encoded by upregulated DEGs. Among these, ellipses with black borders denote the top 10 downregulated DEGs in glaucoma. The ellipses with blue border were represented the top 10 upregulated DEGs. DEGs, differentially expressed genes.

TFs of the top 20 DEGs and TF regulation network

A total of 250 TF binding associations were identified, including multiple binding sites of the same TF in a gene. Among these, 36 TFs were involved. In Table II, the top six TFs with the most downstream genes (the top 20 differentially expressed genes) were included: PAX4; solute carrier family 22 member 1 (1-Oct); hepatocyte nuclear factor 4 α (HNF-4); NK2 homeobox 5 (Nkx2-5); PAX6; and ELK1, ETS transcription factor (Elk-1).

Table II. Top six TFs with the highest number of downstream regulatory genes and their target genes.

Table II.

Top six TFs with the highest number of downstream regulatory genes and their target genes.

TF name	Number of regulated genes	Regulated genes
Pax-4	11	HES4, ITM2C, ANKH, HLA-DPA1, RNF122, LCN2, CCDC3, NELL2, MSLN, CHST6, KRT19
1-Oct	11	ITM2C, MASP1, HOPX, ANKH, HLA-DPA1, OSR2, RNF122, CCDC3, HBB, NELL2, MSLN
HNF-4	11	ITM2C, MASP1, HLA-DPA1, LCN2, MAOA, HBB, FADS2, NELL2, MSLN, CHST6, KRT19
Nkx2-5	9	PPT1, HOPX, RNF122, CCDC3, FADS2, NELL2, MSLN, CHST6, KRT19
Pax-6	8	ITM2C, MASP1, ANKH, OSR2, CCDC3, NELL2, GCNT3, KRT19
Elk-1	5	ITM2C, MASP1, ANKH, LCN2, CCDC3, FADS2, CHST6

[i] TF, transcription factor.

Cytoscape software was used to identify the TFs and the differentially expressed regulatory network of their target genes (Fig. 5). There were 61 nodes, including 44 transcription factors and 20 differentially expressed target genes and 132 edges identified. The eight target genes with the highest degrees were: Coiled-coil domain containing 3 (degree=14); neural tissue-specific epidermal growth factor-like repeat domain-containing protein (NELL2; degree=13); hes family bHLH transcription factor 4 (degree=10); mesothelin (degree=10); integral membrane protein 2C (degree=9); fatty acid desaturase 2 (degree=9); mannan binding lectin serine peptidase 1 (degree=8); and hemoglobin subunit b (HBB; degree=8).

Figure 5. Glaucoma-specific transcription factors regulation network. Blue diamonds represent transcription factors, and ovals represent the top 20 differentially expressed genes. Red ovals denote upregulated genes and green ovals represent downregulated genes.

Validation of DEGs in GEO GSE9944 dataset

To define the key genes with key roles in glaucoma, the GEO GSE9944 repository was searched for high throughput gene expression data and hybridization arrays, chips and microarrays. As demonstrated in Fig. 6, the genes CREB1 and FN1 were upregulated, and the genes KRT19, LCN2 and PAX6 were downregulated. These results were consistent with the integrated analysis in the GSE27276 and GSE4316 datasets.

Figure 6. Validation of the expression levels of different genes in glaucoma compared with normal control based on the Gene Expression Omnibus GSE9944 database. (A) CREB1, (B) FN1, (C) KRT19, (D) LCN2 and (E) PAX6. The expression of these selected genes was consistent with the integrated analysis in the GSE27276 and GSE4316 datasets. CREB1, cAMP responsive element binding protein 1; FN1, fibronectin 1; KRT19, keratin 19; LCN2, lipocalin 2; PAX6, paired box 6.

Discussion

Glaucoma is the leading cause of irreversible blindness worldwide and there is no effective treatment at present (27). In the present study, integrated analysis was performed using data obtained from the GEO database. KEGG, GO and other biological information databases, and R analysis tools were used to analyze the DEGs. A total of 1,935 DEGs in glaucoma (951 genes were upregulated, 984 genes were downregulated) were obtained. Critical signaling pathways that affected the pathogenesis of glaucoma, including ‘eye development’ (FDR=0.00415533) and ‘visual perception’ (FDR=0.00713283) were identified. In addition, based on the promoter sequence of DEGs obtained from UCSC, a TF regulatory network was constructed using the match tool of the TRANSFAC website to obtain the corresponding TFs.

LCN2, encoded by lcn2 gene, is a neutrophil gelatinase-associated lipocalin. The 25-kD LCN2 protein is able to bind small lipophilic substances, including bacteria-derived lipopolysaccharide and formylpeptides (28). The protein additionally functions as a modulator of inflammation. Khalyfa et al (29) suggested that LCN2 was highly upregulated in glaucoma. In the results of the present study, LCN2 was downregulated and the validation result in the GSE9944 dataset was consistent with the analysis. A previous study additionally demonstrated that the coadministration of a reversible monoamine oxidase A (MAOA) inhibitor with epinephrine may be useful for patients with glaucoma (28). In the results of the present study, MAOA was downregulated.

CREB1 is a transcription regulator activated in response to harmful stress stimuli, including hypoxia and oxidative stress, and is involved in the cellular defense against these stresses (30). High IOP, optic nerve damage and visual field defects are the primary pathological features of glaucoma, with high IOP being the most common (27). CREB has been demonstrated to promote neuronal survival (31), and a previous study indicated that CREB has a neuroprotective effect against hydrogen peroxide-induced retinal ganglion cell death via 2 downstream cell survival genes, including brain-derived neurotrophic factor and apoptosis regulator Bcl-2 (32). In a study by Yasuda et al (12), an in silico pathway analysis conducted in rats with axonal injury induced by transecting the right optic nerve, CREB1 was the most significant upstream regulator. In the present study, CREB1 (degree=11) was the in the top 3 genes with the highest degrees in the PPI network. In the KEGG analysis, CREB1 was enriched in the function of Huntington's disease.

FN1 is a glycoprotein within the extracellular matrix, which has a crucial role in cell adhesion, growth, migration and differentiation (33,34). Cellular FN is secreted in a soluble form by a variety of different cell types, including fibroblasts, chondrocytes, macrophages and epithelial cells, and is later assembled into an insoluble extracellular matrix (35). FN has a multimodular structure organized into functional domains that interact with multiple binding partners, including integrins, heparin sulfate proteoglycans, collagen, glycosaminoglycans, proteoglycans, heparin, fibrin and bacteria (36). A previous study by Anshu et al (2) identified that FN1 (P=0.0003) was expressed at a significantly increased level in aqueous humor samples compared with controls. In the present study, FN1 was the top protein with the highest degree in the PPI network. Additionally, in the KEGG analysis, FN1 was enriched in the function of ‘Focal adhesion’.

NELL2 is a secreted glycoprotein that is predominantly expressed in neural tissues. NELL2 is additionally involved in promoting the neuronal survival required for the formation of a sexually dimorphic nucleus of the preoptic area in male rats (37). In the TF regulation network in the present study, NELL2 (degree=13) was the in top 2 genes and was regulated by PAX4, 1-Oct, HNF-4, Nkx2-5 and PAX6. Additionally, NELL2 was among the top 40 DEGs. These results indicated that NELL2 may have a key role in the pathogenesis of glaucoma.

The OMIM database was also utilized in the present study, and it was identified that PAX6 was key gene. PAX6 is a member of the paired box gene family, which encodes a transcriptional regulator involved in oculogenesis and other developmental processes. Halder et al (38) suggested that the ectopic expression of Drosophila PAX6 induces ectopic eye development. Wawersik and Maas (39) reviewed the role of PAX6 and other genes in vertebrate and fly oculogenesis; the study demonstrated that regulators of eye development are conserved across large evolutionary distances. The expression patterns and conserved functional domains suggested that Pax6 may have important roles in vertebrate eye formation. Fujimura et al (40) demonstrated that PAX6 appears to function as a pleiotropic regulator, directing the development of ocular tissues in concert with the signaling pathway and, concomitantly, regulating the expression of structural components of the eye, including shielding pigment. In the present study, PAX6 was downregulated in the patients with glaucoma. In addition, in the PPI network, PAX6 was one of the top eight genes with the highest degrees.

In conclusion, six DEGs (LCN2, MAOA, HBB, PAX6, FN1 and CREB1) were identified to be involved in the process of glaucoma. From the two GEO datasets analyzed, 1,935 DEGs (951 upregulated and 984 downregulated genes) were identified between glaucoma and normal controls. DEG validation in the GEO GSE9944 dataset was consistent with the integrated analysis of the GSE27276 and GSE4316 datasets. These results may contribute to the elucidation of novel potential biomarkers, reveal the underlying pathogenesis and identify novel therapeutic targets for the treatment of glaucoma.

Acknowledgements

Not applicable.

Funding

The present study was supported by Science and Technology Development Plan of Jining [grant no. (2016) 56–31] and Project Medical and Health Technology Development Project of Shandong (grant no. 2017WS528).

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Authors' contributions

JF and JX conceived the study and wrote the manuscript. JX performed the data analyses. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References