Skip to main content
Download PDF

Gene network analyses unveil possible molecular basis underlying drug-induced glaucoma

BMC Medical Genomics volume 14, Article number: 109 (2021) Cite this article

Abstract

Background

Drug-induced glaucoma (DIG) is a kind of serious adverse drug reaction that can cause irreversible blindness. Up-to-date, the molecular mechanism of DIG largely remains unclear yet due to the medical complexity of glaucoma onset.

Methods

In this study, we conducted data mining of tremendous historical adverse drug events and genome-wide drug-regulated gene signatures to identify glaucoma-associated drugs. Upon these drugs, we carried out serial network analyses, including the weighted gene co-expression network analysis (WGCNA), to illustrate the gene interaction network underlying DIG. Furthermore, we applied pathogenic risk assessment to discover potential biomarker genes for DIG.

Results

As the results, we discovered 13 highly glaucoma-associated drugs, a glaucoma-related gene network, and 55 glaucoma-susceptible genes. These genes likely played central roles in triggering DIGs via an integrative mechanism of phototransduction dysfunction, intracellular calcium homeostasis disruption, and retinal ganglion cell death. Further pathogenic risk analysis manifested that a panel of nine genes, particularly OTOF gene, could serve as potential biomarkers for early-onset DIG prognosis.

Conclusions

This study elucidates the possible molecular basis underlying DIGs systematically for the first time. It also provides prognosis clues for early-onset glaucoma and thus assists in designing better therapeutic regimens.

Peer Review reports

Background

Glaucoma is a severe disease that affects more than 70 million people worldwide. As a leading cause of irreversible blindness, glaucoma’s clinical manifestations include optic disc atrophy, visual field defect, and visual acuity loss; it brings great inconvenience and desperate darkness to the patients [1]. Drug-induced glaucoma (DIG) is a special kind of glaucoma, which may be caused by many commonly used drugs, such as dexamethasone, paclitaxel, fluoxetine, and perphenazine in normal chemotherapy [2, 3]. Generally, the DIGs can be classified into two types: the open-angle glaucoma (OAG) with elevated intraocular pressure (IOP) which is mainly caused by the steroids, such as betamethasone, prednisolone, and dexamethasone; and the closed-angle glaucoma (CAG) with blocked trabecular network outflow path which is mainly caused by the non-steroidal drugs such as anticholinergic drugs and epinephrine drugs. Some antitumor drugs, such as docetaxel and paclitaxel, can also induce OAG via unclear internal mechanisms [2,3,4]. As of December 2019, a total of 13,961 cases (including 12,195 serious cases) of DIGs have been reported in the FDA Adverse Event Reporting System (FAERS) of the United States.

In recent years, extensive attention has been paid to discovering the molecular mechanisms beneath glaucoma. Both the animal model and the Genome-Wide Association Studies (GWAS) were adopted to explore susceptible pathways or genes beneath glaucoma and obtained some encouraging achievements. For instance, Weinreb linked glaucoma with multiple pathways of lipid metabolism, inflammation, autoimmunity, and oxidative stress [1]. Razeghinejad attributed forward movement of the iris–lens diaphragm to the closed-angle glaucoma [3]. Clark proposed glucocorticoid-induced OAGs likely led by inhibition of matrix metalloproteinases [5]. Danford suggested the primary OAGs related to MYOC, OPTN, CAV1, and CAV2 [6]. Regretfully, these efforts failed to identify a strong genetic factor in glaucoma development. The surgery-simulated glaucoma in the animal model experiment is to some extent different from the real clinical case, which may bring unpredictable errors. Taking the advantages of high-throughput technology, the GWAS studies discovered a list of genes and genomic mutations susceptible to DIG on a large scale [7]. We summarized these efforts in Additional file 1: Table S1, in which information of four potential DIG-associated genes (MYOC, GPNMB, HCG22, and TSPO) were given in brief. Regretfully, the GWAS applications are largely limited by sample size and financial budgets. Moreover, many of the susceptible genes identified by the GWAS studies have relatively small penetration to glaucoma because of the inherent limitations of GWAS itself, for instance, the inability to fully explain the genetic risk of common diseases and the difficulty in figuring out the true causal associations [8]. Therefore, more susceptible genes are awaiting discovery from the larger DIG cohorts, and at the same time, well-defined phenotypes are expected to delineate the complex genetic architecture of glaucoma [9].

In the present study, we attempt to integrate and analyze the heterogeneous data of clinical records and transcriptomes to identify the drugs that are highly associated the DIGs, unveil the gene interactions underlying glaucoma development, extract the pivot pathways and genes underlying glaucoma, and eventually discover the potential gene biomarkers for DIG prognosis.

Methods

Data source and preprocessing

The drug-ADR relations

The drug-ADR relations were derived from the Adverse Drug Reaction Classification System (ADReCSv2.0) [10]. Overall, 173 distinct glaucoma-inducing drugs were extracted from the ADReCS, as well as their comprehensive drug information and ADR information. The drug indication information was subject to the Anatomical Therapeutic Chemical Classification System (ATC Code).

The drug-treated gene expression profiles

The drug-treated gene expression profiles were downloaded from the LINCS project [11]. Overall, 1,319,138 gene expression profiles were obtained from the LINCS database, covering 25,149 chemicals, 12,328 genes, 76 cell types, and 5,420 combinational conditions (cell line, drug dosage, and exposure time). For every drug-treated profile, the gene expression perturbation to the control was determined by the LINCS L1000 project using the Characteristic Directions (CD) geometrical approach [12] and as well the significance (p value) of perturbation was calculated [11]. In addition, the top perturbed genes were elected as the differentially expressed genes (DEGs). The profiles were then preprocessed as following: (1) excluded the experiments without biological replicates; (2) for every drugs, chose the mostly perturbed profiles (with minimum p value and p < 0.01) as the representatives. After the data preprocessing, overall 1114 gene expression profiles were retained for later analysis, covering 299 chemicals (272 approved, 4 experimental, and investigational drugs), 12,328 genes, 9 cell types, and 27 combinational experiment conditions (drug dosage and exposure time) (Additional file 2: Table S2).

Mining drug-glaucoma association

To obtain reliable drug-glaucoma associations, we downloaded 9,772,360 adverse drug event (ADE) reports, dating from Jan 2004 to Dec 2018, from the FDA Adverse Event Reporting System (FAERS) of the United States. Excluding the reports from unreliable sources and non-single active ingredient drugs, a total of 2,069,653 reports were used for association analysis, including 1095 distinct single-ingredient drugs and 10,287 adverse reactions. The ADR terms were standardized by referring to the ADReCS. Subsequently, we calculated the odds ratio (OR) for every drug-ADR pairs following the standard procedure. This analysis involved 106 distinct drug-glaucoma associations after ADR standardization.

Tissue propensity analysis

We manually mapped the nine studying cell types against the System Organ Class (SOC) of the ATC system, each cell type corresponded to only one organ. As the result, four organs were involved (Additional file 3: Table S3). Of these organs, the ATC code of a drug at the SOC level were taken as its expected site-of-action (SOA). For every drugs involved in this analysis, we calculated the average expression changes of its DEGs in response to the drug treatment. The average perturbations (the average expression changes of DEGs) of drugs in SOA and non-SOA groups in each cell type were statistically analyzed.

The weighted gene co-expression network analysis (WGCNA)

The weighted gene co-expression network analysis (WGCNA) was conducted on 13 selected drug-treated gene expression profiles using the R package “WGCNA” [13], and the soft threshold for network construction was set to 5 (with scale-free coefficient R2 = 0.88). Upon the gene network, the eigengene dendrogram was constructed using the average linkage hierarchical clustering and the dynamic tree-cutting algorithm [14]. A dissimilarity threshold of 0.3 (or similarity threshold of 0.7) was set for the eigengene dendrogram to reduce network complexity by merging similar modules. The genes in the same module were taken as the co-expressed genes. For each gene module, the module eigengenes (MEs) were the first principal component of the expression matrix.

Identification of glaucoma-associated genes

For the selected gene module, we quantitatively measured the module membership (MM) between MEs by the Pearson Correlation analysis. The highly connected genes (MM > 0.95) were selected as core genes for subsequent protein–protein interaction (PPI) analysis. The PPI data for highly connected genes were obtained from the STRING database (v10.5, species: Homo sapiens) by satisfying PPI enrichment significance p < 0.01 [15]. Upon the PPI data, we re-constructed the PPI network using the Cytoscape software (version: 3.7.1) [16], from which we extracted the hub genes (Maximal Clique Centrality ≥ 5 and Betweenness centrality > 0) as the potential glaucoma-associated genes via the CytoHubba plugin.

Functional enrichment analysis

Gene Ontology (GO) and KEGG pathways enrichment analysis were conducted on the genes of interest using the R package ClusterProfiler [17]. Fisher’s exact test was used to evaluate the significance of the biology processes or pathways. The false discovery rate (FDR) correction method (Benjamini and Hochberg) was used to evaluate the multiple tests in the enrichment analysis.

The permutation test

We performed the permutation test to measure the significance of gene expression differences under the treatment of 13 glaucoma-associated drugs and 63 non-glaucoma-inducing drugs. The drugs of both drug groups all belonged to the ATC category "S" (Sensory organs). For each of the studying genes, the permutation test was demonstrated independently on a null hypothesis (i.e. no gene expression difference between two drug groups). We simulated the drug group composition by resampling the group members 10,000 times randomly, calculated the gene expression difference between drug groups, and determined the distribution of gene expression difference. At the same time, we calculated the real gene expression difference between two drug groups for each gene. Subsequently, we located the real difference at the random expression difference distribution, upon which the significance p value was determined by counting the proportion of samples less than the real difference value. A p < 0.05 will reject the null hypothesis, indicating there exists a significant gene expression difference between the two drug groups.

Pathogenic risk assessment

The pathogenic risk assessment was demonstrated in this study by calculating the odds ratio (OR) value of the selected gene to glaucoma. The odds ratio analysis took the 77 gene expression profiles as the independent variables x and the binary consequence of glaucoma (1 for glaucoma and 0 for non-glaucoma) as the response variable Y [12]. We assumed there existed a linear relationship between the independent variables x and the log-odds of \({\text{Y}} = 1\), denoted as \(p = {\text{P}}\left( {{\text{Y}} = 1} \right)\). This linear relationship can be determined by:

$$\ell = \log_{e} \frac{p}{1 - p} =\upbeta _{0} +\upbeta _{i} x_{i}$$
(1)

where ℓ stands for the log-odds, \({\upbeta }_{i}\) stands for the parameter of the Logistic Regression model. According to the relationship of the response variable Y and independent variable x, we drew the receiver operating characteristic curve (ROC curve) and calculated the area under the curve (AUC) for each variable x. In this study, the logistic regression was performed using the R function glm (R version v3.60), and the AUC was calculated by R function colAUC from caTools (v1.18.0).

Results

Identify the glaucoma-associated drugs

Overall, 173 glaucoma-inducing drugs were derived from the ADReCSv2.0 database. These drugs can induce glaucoma at different occurrences in patients as estimated in the drug label. To determine the drug-glaucoma associations quantitatively, we analyzed 9,772,360 real-world historical ADEs of the FAERS using the odds ratio analysis and the Fisher test. A list of 25 significant (OR > 2 and p < 0.05) drug-glaucoma associations were summarized in Additional file 4: Table S4. In the list, 13 drugs exhibited strong associations with glaucoma (ORs > 10 and p < 0.01). Additional analysis found that drugs of the ATC code "S" (Sensory organs) category were more significantly associated with glaucoma compared to other non-S category drugs (p < 0.02) (Fig. 1). The average ORs for S category drugs and the non-S category drugs were 21.13 and 14.74, respectively.

Fig. 1
figure 1

Identification of significant drug-glaucoma associations. a The association strength (odds ratio) of 71 glaucoma-inducing drugs among 10 ATC categories. b The association strength comparison of 10 drugs in ATC category “S” (Sensory organ) and 61 drugs in other categories. The drug-glaucoma association strength were determined upon 2,069,653 historical ADEs of the FARES database using the odds ratio analysis. One-sided Wilcox test was used to determine the significance of the difference between these two groups. The star symbol (*) stands for p value < 0.05. The drugs of category S had higher association strength than that of other drug categories

Full size image

Furthermore, we conducted a comparative analysis of genome-wide gene expression changes over different cell lines in response to drug treatment. The result manifested that the gene expressions were more perturbed in the constituted cells of the site-of-actions (SOAs) of drugs than in the cell lines of other tissues/organs (Fig. 2). This suggested that many high occurrence ADRs were likely caused by the augment effects of gene expression perturbation. Therefore, we speculated that DIGs might tend to occur on the sensory organs, particularly the eyes. Accordingly, 13 drugs were selected for later mechanistic analyses by satisfying both criteria of belonging to the S category and having a strong association with glaucoma (Additional file 5: Table S5).

Fig. 2
figure 2

The tissue propensity of drug-induced glaucoma. The X-axis stands for the nine different cell types and Y-axis stands for average perturbation of the DEGs. The pink box stands for the expected Site-of-Action (SOA) of the drugs. The blue box stands for the tissues/organs other than the SOA of the drugs. One-sided Wilcox test was used to determine the significance of the difference between SOAs and other tissues/organs. The star symbol (*) stands for p < 0.05. The results manifested that the drugs intended to interfere gene expression mostly at its SOA

Full size image

Before moving forward, we made a structural similarity measurement of the selected 13 glaucoma-associated drugs using the Tanimoto coefficient, in which the chemical structures were represented by the binary fingerprint value. The measurement revealed a low similarity between drugs (the average Tanimoto Coefficient = 0.43). As well, we evaluated the dispersion degree of physicochemical properties for the 13 drugs using the coefficient of variation (CV). The results also supported the high variability of drugs (average CV = 0.358). Both analyses manifested that the 13 glaucoma-associated drugs were highly diverse in structure and physicochemical properties; and as a consequence, it is unlikely for the drugs to induce glaucoma in a simple way of targeting the same protein.

Unveil the gene interactions underlying the glaucoma development

To unveil gene interactions underlying the DIGs, we carried out WGCNA analyses on gene expression profiles under the treatment of the 13 selected glaucoma-associated drugs as described in “Methods” section. The analysis extracted 72 gene modules (subnetworks) from the whole gene co-expression network via average linkage hierarchical clustering and dynamic tree-cutting algorithm. To reduce the network complexity, we set a dissimilarity threshold of 0.3 (or similarity value = 0.7) to merge gene modules in the dendrogram of module eigengenes (MEs) and eventually obtained 27 modules (Additional file 7: Figure S1). These 27 gene modules or subnetworks represented the functional gene groups in response to the treatment of 13 glaucoma-associated drugs. Subsequent function enrichment analyses on these 27 gene modules identified one module (involving 5,968 genes) potentially linked with DIG. This gene module contained the biological function terms related to sensory organ such as "visual perception" (GO:0007601), "sensory perception of light stimulus" (GO:0050953), and "Phototransduction" (hsa04744). We further compared the functional enrichment results to two prior studies of glaucoma [6, 18]. As the results, six GO terms (in top 15) and 14 KEGG pathways (in top 20) were found mutual to both prior studies and this study, which affirmed that the selected module was the potential DIG-related gene module (Additional file 8: Figure S2).

Identify glaucoma-associated genes

To identify potential DIG-associated genes/proteins in the DIG-related gene module, we carried out serial network analyses step by step. Firstly, we calculated the intra-module network connectivity for all genes in the DIG-related module. The calculation determined 326 highly connected genes by satisfying the criteria of connectivity > 0.95. These highly connected genes likely played centralized roles in the DIG-related gene interaction network. Subsequently, we searched these 326 genes against the STRING database, which identified 495 additional gene/protein interactions by meeting the criteria of the PPI enrichment significance p < 0.01. Upon the 326 genes and corresponding 495 interactions, we constructed a PPI network using the Cytoscape software. On basis of the PPI network, we further extracted a core subnetwork of the most enrichment (having the most connectivity) using the Cytoscape tool CytoHubba. This core subnetwork included 55 hub genes, which were most likely associated with DIG. Of them, 13 genes were previously reported to be associated with the sensory functions (Fig. 3). These 13 glaucoma-associated genes can be divided into three functional groups: phototransduction-related genes (GNGT1, OPN1SW, PDE6A, PDC, CRX, CNGB1, GUCY2F), glutamate receptor (GR) genes (GRIN2B, GRIK3, GRIK4), and 5-hydroxytryptamine receptor (HTR) genes (HTR1A, HTR1E, HTR7). Their connections with sensory functions and glaucoma were concluded in Table 1. It is a reasonable assumption that the 55 glaucoma-associated genes may react more actively to the treatment of the 13 glaucoma-associated drugs than that of the 63 non-glaucoma-inducing drugs. To prove the assumption, we carried out the permutation test. As the results, all of the 55 genes showed lower, 50 significantly (p < 0.05) and 5 insignificantly, gene expressions in response to the glaucoma-associated drugs. Furthermore, we observed a positive correlation (Pearson correlation coefficient = 0.56, p = 1.04e−5) between the differential significance (−log10p_value) and glaucoma-gene association strength (OR) of these 55 genes (Additional file 9: Figure S3). For the 13 reported glaucoma-associated genes, 12 genes had significantly (p < 0.05) lower expressions under the treatment of the glaucoma-associated drugs than that of non-glaucoma-associated drugs. The remaining one PDE6A showed a clearly lower expression but statistically insignificant (p = 0.061) (Fig. 4). All these results further supported that the down-regulation of the 55 genes likely answered for the DIGs.

Fig. 3
figure 3

The glaucoma-related gene interaction network. This network was constructed in basis of 326 highly connected genes and corresponding 495 protein–protein interactions using the CytoScape tool. The purple nodes in the middle circle were 42 glaucoma-associated genes, and the innermost nodes were 13 DIG-associated genes supported by literature surveillance

Full size image
Table 1 The information of 13 glaucoma-associated genes
Full size table
Fig. 4
figure 4

The permutation test of 13 glaucoma-associated genes. The test was conducted on 13 glaucoma-associated drugs and 63 non-glaucoma-inducing drugs. The density plot illustrated the distribution of gene expression difference after 10,000 resamplings. The red line stands for the position of real expression difference on the distribution curve. The boxplot illustrated the average gene expression changes under the treatment of two drug groups. The symbol * stands for the significance p value < 0.05 and ** stands for p value < 0.01

Full size image

The down-regulation of these genes by glaucoma-associated drugs were also reported in previous studies. For instance, dexamethasone and triamcinolone are commonly used in clinical practice to treat a wide range of retinal pathologies and both can induce glaucoma [26]. In mice, OPN1SW, GRIN2B, HTR1A, and HTR7 were found downregulated after retina postintravitreal injections of triamcinolone [27]. CNGB1, HTR1A, HTR7, OPN1SW, and CRX were downregulated by dexamethasone in the retinal pigment epithelium of mouse [28].

Discover potential biomarker genes for DIG

A good gene biomarker should be able to indicate disease incidence and as well be highly sensitive to the right drug treatment. The above analyses have confirmed the association of 55 genes with glaucoma. Here, we further quantitatively examined the pathogenic risk of these genes to glaucoma and their sensitivity to the glaucoma-associated drugs. For this purpose, we calculated the odds ratio (OR) for every 55 genes to glaucoma. And meanwhile, we determined the AUC for every gene alone in differentiating glaucoma-associated drugs and non-glaucoma-inducing drugs. As illustrated in Fig. 5 and Additional file 6: Table S6, 9 out of 55 genes exhibited highly risky to glaucoma and highly sensitive to glaucoma-associated drugs (OR > 6 and AUC > 0.7). These nine genes (OTOF, CRX, GLRA1, XCR1, TAS2R13, PDC, GIPR, GNAT3, and TACR1) could constitute a gene panel for diagnosing DIG. In particular, the otoferlin gene (OTOF) exhibited the riskiest to glaucoma with OR = 9.727, and it also showed remarkable ability in discriminating the glaucoma-inducing drugs from the non-glaucoma-inducing drugs with an AUC = 0.755. Hence, OTOF may serve as a potential biomarker gene for indicating DIGs.

Fig. 5
figure 5

The pathogenic risk assessment of 55 glaucoma-associated genes. The bar length is positively proportional to the odds ratio (OR) value of the gene-to-glaucoma; the bar color is proportional to the AUC value for the gene ability in differentiating the glaucoma-associated drugs and non-glaucoma-inducing drugs. The color of gene symbol corresponds to different gene group: green stands for the phototransduction-related gene, yellow stands for the Glutamate receptor gene, red stands for the 5-hydroxytryptamine receptor gene, and black stands for the others

Full size image

Discussion

The drug-induced glaucoma (DIG) has been encountered in normal chemotherapy for a long time but has not yet been fully understood. One of the major reasons is that glaucoma-induced drugs cannot provide suggestive chemical or protein target clues. In this study, we identified 13 strong glaucoma-associated drugs via mining 9,772,360 historical ADEs in the FAERS. With no surprise, these 13 drugs did not show significant structural or physicochemical similarity, which made it hard to explain the DIGs in a conventional way of drug-target interactions. Alternatively, we analyzed the mutual gene expression patterns treated by the 13 DIG-associated drugs using the WGCNA method, via which we uncovered a potential DIG-related gene module (network) and 55 potential glaucoma-associated genes. Subsequent network analyses and expression difference analysis further consolidated that the 55 genes might play central roles in the glaucoma development. Of these glaucoma-associated genes, 13 genes were literature-reported and can be divided into three functional gene groups according to their possible roles in glaucoma development, including seven phototransduction-related genes, three GR genes, and three HTR genes. According to the functional gene groups, we proposed here three possible routes, but not limited to, to induce glaucoma: phototransduction dysfunction, intracellular calcium homeostasis disruption, and retinal ganglion cell (RGC) death.

In clinical practices, glaucoma was always accompanied with several symptoms, such as elevated intraocular pressure (IOP), visual field defect, and optic nerve damage [1, 3]. Phototransduction is one of the main processes in vision formation. Of the seven phototransduction-related genes, GNGT1 (the guanine nucleotide-binding protein G (T) subunit gamma-T1) is specifically expressed in the retinal rod outer segment cells, participating in the signal transduction process of cyclic GTP-specific phosphodiesterase as a modulator or transducer of light transduction. Down-regulation of GNGT1 may result in blockage of light conduction, leading to visual field defect and vision loss. A previous study has shown that elimination of GNGT1 could cause a dramatic decrease of the detected light signal in the intact mouse rods, a striking decline in the rod visual sensitivity, and eventually severe impairment of the nocturnal vision [29]. For other phototransduction-related genes CNGB1, CRX, GUCY2F, OPN1SW, PDC, and PDE6A, they may participate in the regulation of the state transition and the phototransduction cascade via interacting with GNGT1 closely.

Glutamate receptors (GRs) and 5-hydroxytryptamine receptors (HTRs) are pivot proteins in maintaining the homeostasis of various neurotransmitters, the transmission of visual information, and the activity of nerve cells. A prior study found that abnormal increase of glutamate concentration in the vitreous bodies of humans and monkeys with glaucoma [30]. Over-activation of NMDA glutamate receptors (GRIN2B) affected the intracellular calcium homeostasis and induced cell apoptosis via caspases and PARP [31, 32], particularly the death of RGCs. The RGC was thought to be the most susceptible cell type to clinical glaucoma [33], and the RGC apoptosis was often observed in human and experimental animal models of glaucoma [34]. For the kainic acid (KA) glutamate receptors (GRIK3 and GRIK4), they were acknowledged to mediate glutamate-mediated neural signals between photoreceptors and bipolar cells [24, 35]. Abnormal regulation of GRIKs may lead to the defect transmission of visual information. For the three HTRs (HTR1A, HTR1E, and HTR7), limited works have been undertaken to link them with glaucoma. However, prior work used HTR2 (serotonin-2 receptor) agonists to reduce ocular hypotension in the chronic glaucoma model [36], which hint the HTRs might be involved in the IOP and glaucoma.

For the 55 glaucoma-associated genes, we further measured their pathogenic risks to glaucoma, as well as their sensitivity to the glaucoma-associated drugs. These intentions quantitatively ranked gene susceptibility to glaucoma, which finally suggested a nine-gene panel (OTOF, CRX, GLRA1, XCR1, TAS2R13, PDC, GIPR, GNAT3, and TACR1), which showed remarkable capability in differentiating glaucoma onset or not. OTOF encodes an integral membrane protein, which is implicated in a late stage of exocytosis and plays a ubiquitous role in synaptic vesicle trafficking. It was known to be a cause of neurosensory nonsyndromic recessive deafness via modulating γ-aminobutyric acid (GABA) activity, the metabolite of an excitatory neurotransmitter glutamate [37]. However, it is suspected that defective OTOF activity would markedly interfere GABA and glutamate metabolism [37], for example, causing abnormal accumulation of glutamate in the retina. The excessive glutamates would likely induce glutamate excitotoxicity to retinal neurons by overstimulation of glutamate receptors [38, 39]. Both the cone-rod homeobox (CRX) gene and the phosducin (PDC) gene are phototransduction-related genes. CRX is a photoreceptor-specific transcription factor that is necessary for the maintenance of normal cone and rod function. PDC encodes a phosphoprotein in the rod cells of retina, which may participate in the regulation of visual phototransduction or in the integration of photoreceptor metabolism. Several previous studies indicated that defects of CRX or PDC might answer for photoreceptor degeneration, leading to cone swelling and loss of photoreceptor cells, and eventually causing blindness [40,41,42,43]. The gastric inhibitory polypeptide receptor (GIPR) can stimulate insulin release in the presence of elevated glucose; variants of GIPR were found linked with diabetes-related primary OAG [44]. For the remaining five biomarker genes, no substantial evidences have been raised till now to link with glaucoma directly yet, which desires further experimental investigation.

This work has its limitations. First of all, the gene network study takes the advantage by bypassing the conventional challenge in discovering drug-target interactions beneath the DIGs; at the meanwhile, it ignores the diverse upstream mechanisms to trigger the DIGs, for instance, the genetic characteristics, the pathogenic conditions, and the drug types such as corticosteroids and anticholinergic drugs. Moreover, the results of this work could be limited by the data availability that some drugs were excluded in the WCGNA analysis due to the lack of sufficient drug-treated gene expression profiles. Nevertheless, this work provides a new insight into the systematic understanding of drug-induced glaucoma from a gene interaction perspective. The same strategy can be easily applied to mechanistically investigate other severe adverse drug reactions. It will also help prevent drug-induced glaucoma in clinical practice by suggesting the potential biomarkers for glaucoma prognosis.

Conclusions

In summary, the development of DIG is a sophisticated process that involves multiple genes and pathways. Current molecular understanding of DIGs is mostly focused on sporadic genes or pathways, which show undetermined impacts on glaucoma. In this work, we develop a gene network analysis strategy to illustrate three possible molecular routes to induce glaucoma. These three routes are not isolated but connected in a glaucoma-related gene network. For the first time, we propose a list of genes highly susceptible to glaucoma. In particular, OTOF shows promising prognostic potential to drug-induced glaucoma in clinics.

Availability of data and materials

The drug-ADR relations were collected from the Adverse Drug Reaction Classification System (http://bio-add.org/ADReCS/download.jsp). The drug-treated gene expression profiles were downloaded from the GEO public database (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE70138). The adverse drug event (ADE) reports were downloaded from openFDA (https://open.fda.gov/data/downloads/).

Abbreviations

DIG:

Drug-induced glaucoma

WGCNA:

Weighted gene co-expression network analysis

OAG:

Open-angle glaucoma

IOP:

Intraocular pressure

CAG:

Closed-angle glaucoma

FAERS:

FDA Adverse Event Reporting System

GWAS:

Genome-Wide Association Studies

CD:

Characteristic Direction

DEG:

Differentially expressed gene

ADE:

Adverse drug event

OR:

Odds ratio

ME:

Module eigengene

MM:

Module membership

PPI:

Protein–protein interaction

GO:

Gene Ontology

FDR:

False discovery rate

ROC:

Curve receiver operating characteristic curve

AUC:

Area under the curve

CV:

Coefficient of variation

HTR:

5-Hydroxytryptamine receptor

RGC:

Retinal ganglion cell

GR:

Glutamate receptor

References

  1. Weinreb RN, Aung T, Medeiros FA. The pathophysiology and treatment of glaucoma: a review. JAMA. 2014;311(18):1901–11.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  2. Razeghinejad MR, Pro MJ, Katz LJ. Non-steroidal drug-induced glaucoma. Eye (Lond). 2011;25(8):971–80.

    Article  CAS  Google Scholar 

  3. Razeghinejad MR, Myers JS, Katz LJ. Iatrogenic glaucoma secondary to medications. Am J Med. 2011;124(1):20–5.

    Article  CAS  PubMed  Google Scholar 

  4. Boonyaleephan S. Drug-induced secondary glaucoma. J Med Assoc Thai. 2010;93(Suppl 2):S118-122.

    PubMed  Google Scholar 

  5. Clark AF, Wordinger RJ. The role of steroids in outflow resistance. Exp Eye Res. 2009;88(4):752–9.

    Article  CAS  PubMed  Google Scholar 

  6. Danford ID, Verkuil LD, Choi DJ, Collins DW, Gudiseva HV, Uyhazi KE, Lau MK, et al. Characterizing the “POAGome”: a bioinformatics-driven approach to primary open-angle glaucoma. Prog Retin Eye Res. 2017;58:89–114.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Khawaja AP, Cooke Bailey JN, Wareham NJ, Scott RA, Simcoe M, Igo RP Jr, Song YE, et al. Genome-wide analyses identify 68 new loci associated with intraocular pressure and improve risk prediction for primary open-angle glaucoma. Nat Genet. 2018;50(6):778–82.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  8. Tam V, Patel N, Turcotte M, Bosse Y, Pare G, Meyre D. Benefits and limitations of genome-wide association studies. Nat Rev Genet. 2019;20(8):467–84.

    Article  CAS  PubMed  Google Scholar 

  9. Wiggs JL, Hauser MA, Abdrabou W, Allingham RR, Budenz DL, Delbono E, Friedman DS, et al. The NEIGHBOR consortium primary open-angle glaucoma genome-wide association study: rationale, study design, and clinical variables. J Glaucoma. 2013;22(7):517–25.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Cai MC, Xu Q, Pan YJ, Pan W, Ji N, Li YB, Jin HJ, et al. ADReCS: an ontology database for aiding standardization and hierarchical classification of adverse drug reaction terms. Nucleic Acids Res. 2015;43(Database issue):D907-913.

    Article  CAS  PubMed  Google Scholar 

  11. Duan Q, Reid SP, Clark NR, Wang Z, Fernandez NF, Rouillard AD, Readhead B, et al. L1000CDS(2): LINCS L1000 characteristic direction signatures search engine. NPJ Syst Biol Appl. 2016;2:1–12.

    Article  CAS  Google Scholar 

  12. Clark NR, Hu KS, Feldmann AS, Kou Y, Chen EY, Duan Q, Ma’ayan A. The characteristic direction: a geometrical approach to identify differentially expressed genes. BMC Bioinform. 2014;15:79.

    Article  Google Scholar 

  13. Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinform. 2008;9:559.

    Article  CAS  Google Scholar 

  14. Langfelder P, Zhang B, Horvath S. Defining clusters from a hierarchical cluster tree: the Dynamic Tree Cut package for R. Bioinformatics. 2008;24(5):719–20.

    Article  CAS  PubMed  Google Scholar 

  15. von Mering C, Huynen M, Jaeggi D, Schmidt S, Bork P, Snel B. STRING: a database of predicted functional associations between proteins. Nucleic Acids Res. 2003;31(1):258–61.

    Article  CAS  Google Scholar 

  16. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, Amin N, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  17. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Bailey JN, Yaspan BL, Pasquale LR, Hauser MA, Kang JH, Loomis SJ, Brilliant M, et al. Hypothesis-independent pathway analysis implicates GABA and acetyl-CoA metabolism in primary open-angle glaucoma and normal-pressure glaucoma. Hum Genet. 2014;133(10):1319–30.

    Article  PubMed  CAS  Google Scholar 

  19. Kimura A, Singh D, Wawrousek EF, Kikuchi M, Nakamura M, Shinohara T. Both PCE-1/RX and OTX/CRX interactions are necessary for photoreceptor-specific gene expression. J Biol Chem. 2000;275(2):1152–60.

    Article  CAS  PubMed  Google Scholar 

  20. Kacirova M, Kosek D, Kadek A, Man P, Vecer J, Herman P, Obsilova V, et al. Structural characterization of Phosducin and its complex with the 14-3-3 protein. J Biol Chem. 2015;290(26):16246–60.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Sarfare S, McKeown AS, Messinger J, Rubin G, Wei H, Kraft TW, Pittler SJ. Overexpression of rod photoreceptor glutamic acid rich protein 2 (GARP2) increases gain and slows recovery in mouse retina. Cell Commun Signal. 2014;12:67.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  22. Lamb TD. Evolution of phototransduction, vertebrate photoreceptors and retina. Prog Retin Eye Res. 2013;36:52–119.

    Article  CAS  PubMed  Google Scholar 

  23. Mowat FM, Occelli LM, Bartoe JT, Gervais KJ, Bruewer AR, Querubin J, Dinculescu A, et al. Gene therapy in a large animal model of PDE6A-retinitis pigmentosa. Front Neurosci. 2017;11:342.

    Article  PubMed  PubMed Central  Google Scholar 

  24. Borghuis BG, Looger LL, Tomita S, Demb JB. Kainate receptors mediate signaling in both transient and sustained OFF bipolar cell pathways in mouse retina. J Neurosci. 2014;34(18):6128–39.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  25. Zhou X, Zhang R, Zhang S, Wu J, Sun X. Activation of 5-HT1A receptors promotes retinal ganglion cell function by inhibiting the cAMP-PKA pathway to modulate presynaptic GABA release in chronic glaucoma. J Neurosci. 2019;39(8):1484–504.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Phulke S, Kaushik S, Kaur S, Pandav SS. Steroid-induced glaucoma: an avoidable irreversible blindness. J Curr Glaucoma Pract. 2017;11(2):67–72.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Smit-McBride Z, Moisseiev E, Modjtahedi SP, Telander DG, Hjelmeland LM, Morse LS. Comparison of in vivo gene expression profiling of RPE/choroid following intravitreal injection of dexamethasone and triamcinolone acetonide. J Ophthalmol. 2016;2016:9856736.

    Article  PubMed  PubMed Central  CAS  Google Scholar 

  28. Smit-McBride Z, Modjtahedi SP, Cessna CT, Telander DG, Hjelmeland LM, Morse LS. In vivo gene expression profiling of retina postintravitreal injections of dexamethasone and triamcinolone at clinically relevant time points for patient care. Invest Ophthalmol Vis Sci. 2011;52(12):8965–78.

    Article  PubMed  Google Scholar 

  29. Kolesnikov AV, Rikimaru L, Hennig AK, Lukasiewicz PD, Fliesler SJ, Govardovskii VI, Kefalov VJ, et al. G-protein betagamma-complex is crucial for efficient signal amplification in vision. J Neurosci. 2011;31(22):8067–77.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Dreyer EB, Zurakowski D, Schumer RA, Podos SM, Lipton SA. Elevated glutamate levels in the vitreous body of humans and monkeys with glaucoma. Arch Ophthalmol. 1996;114(3):299–305.

    Article  CAS  PubMed  Google Scholar 

  31. Kwong JM, Lam TT. N -methyl-d-aspartate (NMDA) induced apoptosis in adult rabbit retinas. Exp Eye Res. 2000;71(4):437–44.

    Article  CAS  PubMed  Google Scholar 

  32. Lam TT, Abler AS, Kwong JM, Tso MO. N-methyl-d-aspartate (NMDA)–induced apoptosis in rat retina. Invest Ophthalmol Vis Sci. 1999;40(10):2391–7.

    CAS  PubMed  Google Scholar 

  33. Nickells RW, Howell GR, Soto I, John SW. Under pressure: cellular and molecular responses during glaucoma, a common neurodegeneration with axonopathy. Annu Rev Neurosci. 2012;35:153–79.

    Article  CAS  PubMed  Google Scholar 

  34. McKinnon SJ. Glaucoma: ocular Alzheimer’s disease? Front Biosci. 2003;8:s1140-1156.

    Article  CAS  PubMed  Google Scholar 

  35. DeVries SH, Schwartz EA. Kainate receptors mediate synaptic transmission between cones and “Off” bipolar cells in a mammalian retina. Nature. 1999;397(6715):157–60.

    Article  CAS  PubMed  Google Scholar 

  36. Sharif NA. Serotonin-2 receptor agonists as novel ocular hypotensive agents and their cellular and molecular mechanisms of action. Curr Drug Targets. 2010;11(8):978–93.

    Article  CAS  PubMed  Google Scholar 

  37. Wu W, Rahman MN, Guo J, Roy N, Xue L, Cahill CM, Zhang S, et al. Function coupling of otoferlin with GAD65 acts to modulate GABAergic activity. J Mol Cell Biol. 2015;7(2):168–79.

    Article  CAS  PubMed  Google Scholar 

  38. Ishikawa M, Yoshitomi T, Covey DF, Zorumski CF, Izumi Y. Neurosteroids and oxysterols as potential therapeutic agents for glaucoma and Alzheimer’s disease. Neuropsychiatry (London). 2018;8(1):344–59.

    Google Scholar 

  39. Cheung W, Guo L, Cordeiro MF. Neuroprotection in glaucoma: drug-based approaches. Optom Vis Sci. 2008;85(6):406–16.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Assawachananont J, Kim SY, Kaya KD, Fariss R, Roger JE, Swaroop A. Cone-rod homeobox CRX controls presynaptic active zone formation in photoreceptors of mammalian retina. Hum Mol Genet. 2018;27(20):3555–67.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Nork TM, Ver Hoeve JN, Poulsen GL, Nickells RW, Davis MD, Weber AJ, et al. Swelling and loss of photoreceptors in chronic human and experimental glaucomas. Arch Ophthalmol. 2000;118(2):235–45.

    Article  CAS  PubMed  Google Scholar 

  42. Humphries PM, Mansergh FC, Farrar J. Polymorphisms in the phosducin (PDC) gene on chromosome 1q25–32. Am J Hum Genet 1994;55(Suppl.3).

  43. Sokolov M, Yadav RP, Brooks C, Artemyev NO. Chaperones and retinal disorders. Adv Protein Chem Struct Biol. 2019;114:85–117.

    Article  CAS  PubMed  Google Scholar 

  44. Shen L, Walter S, Melles RB, Glymour MM, Jorgenson E. Diabetes pathology and risk of primary open-angle glaucoma: evaluating causal mechanisms by using genetic information. Am J Epidemiol. 2016;183(2):147–55.

    PubMed  Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

This study was supported by a grant from the National Key Research and Developmental Program of China (2018YFC1003601); The funding body had a role in the design of the study and data collection.

Author information

Author notes

  1. Ruo-Fan Ding and Qian Yu have contributed equally to this work and Co-first authors

Authors and Affiliations

  1. State Key Laboratory of Stress Cell Biology, School of Life Sciences, Xiamen University, Xiamen, 361102, Fujian, People’s Republic of China

    Ruo-Fan Ding, Qian Yu, Ke Liu, Juan Du, Hua-Jun Yin & Zhi-Liang Ji

Authors
  1. Ruo-Fan Ding
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Qian Yu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Ke Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Juan Du
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Hua-Jun Yin
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Zhi-Liang Ji
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Z-LJ/KL conceived the study and designed the project. QY/RF-D performed data acquisition and analysis. Z-LJ/R-FD/QY draft the manuscript. DJ/H-JY prepared the final version of the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Zhi-Liang Ji.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors have no conflict of interest to declare.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

. Previously reported DIG-related genes.

Additional file 2

. Detailed information of 1,114 gene expression profiles from the LINCS project.

Additional file 3

. The mapping of the ATC code and cell line.

Additional file 4

. The list of significant drug-glaucoma associations (OR > 2 & p < 0.05).

Additional file 5

. The information of 13 glaucoma-associated drugs.

Additional file 6

. The results of pathogenic risk assessment.

Additional file 7

. The eigengene dendrogram constructed by the WGCNA method.

Additional file 8

. The illustration of gene ontology enrichment.

Additional file 9

. The correlation analysis between the differential significance and glaucoma-gene association strength.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ding, RF., Yu, Q., Liu, K. et al. Gene network analyses unveil possible molecular basis underlying drug-induced glaucoma. BMC Med Genomics 14, 109 (2021). https://doi.org/10.1186/s12920-021-00960-9

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12920-021-00960-9

Keywords

BMC Medical Genomics

ISSN: 1755-8794

Contact us