Identification of Optic Nerve–Related Biomarkers in Primary Open-Angle Glaucoma Based on Comprehensive Bioinformatics and Mendelian Randomization

Abstract

Purpose: Glaucoma is the primary cause of permanent vision loss worldwide. However, the pathogenesis of primary open-angle glaucoma (POAG), the main type of glaucoma, has not yet been completely understood.

Methods: In our study, the POAG cohorts were obtained from the Gene Expression Omnibus (GEO) database (GSE45570). Biomarkers with diagnostic utility for POAG were identified through combining differentially expressed analysis, enrichment analysis, machine learning algorithms, and receiver operating characteristic (ROC) analysis. The regulatory networks (including a competing endogenous RNA (ceRNA) regulatory network and a small molecule compounds-mRNA network) were created. In addition, the Mendelian randomization (MR) analysis was used to identify exposures causally associated with POAG. Finally, the expression of the biomarkers was validated via real-time quantitative polymerase chain reaction (RT-qPCR).

Results: The Gene Ontology (GO) items that the differentially expressed genes (DEGs) between POAG and control groups enriched were relevant to light stimulation and DNA methylation. A total of three light stimulation-related biomarkers (RAB8A, PRG3, and SMAD3) were identified, which had diagnostic value for POAG patients. Besides, the ceRNA regulatory network contained 88 nodes and 93 edges, and a small molecule compounds-mRNA network included 66 nodes and 76 edges. The MR results indicated a causal association between DNA methylation GrimAge acceleration and POAG. Additionally, the results of RT-qPCR revealed that the expression trend of RAB8A was consistent with that of GSE45570.

Conclusions: Taken together, this study provides three light stimulation-related biomarkers (RAB8A, PRG3, and SMAD3) for the diagnosis of POAG, providing scientifically valuable insights for further studies of POAG.

Translational Relevance: Discovering biomarkers that possess diagnostic significance for POAG has the potential to offer new insights into the pathogenesis of POAG and present novel objectives for clinical intervention.

Introduction

Glaucoma is an optic neuropathy disease characterized by visual field defects and the loss of retinal ganglion cells (RGCs), and it is currently recognized as the second most prevalent cause of blindness globally.¹ Based on the anatomical structure of the iris and cornea, glaucoma is categorized as open-angle or closed-angle glaucoma, and, depending on the cause of the disease, it can be categorized as primary or secondary glaucoma, of which primary open-angle glaucoma (POAG) is the most common.² Various factors have been identified as potential risk factors for POAG, including age, family history, race, higher intraocular pressure (IOP), myopia, and systemic disorders.^3,4 Due to the insidious and irreversible nature of POAG, approximately 40% or more of the RGCs have been lost when typical visual field changes appear clinically, resulting in severe visual impairment.⁵ However, to date, the specific mechanism of action leading to optic nerve damage has been unclear for POAG. It is critical to explore the pathogenesis of POAG and develop new diagnostic and therapeutic modalities. 

In recent times, the utilization of microarray and high-throughput technologies has significantly enhanced the efficacy of bioinformatics methods in the identification of differentially expressed genes.^6–9 Consequently, these methods have been widely employed in the investigation of genetic variables associated with diverse ocular illnesses. It has been demonstrated that NEUROD1 serves as a biomarker for POAG retinopathy and plays an essential role in the intricate process of eye development through the utilization of bioinformatics analysis.¹⁰ The identification of eight genes—TOP1, NAV3, RXRB, C1QB, ADAM15, ZNF207, P2RY4, and VAV3—that are involved in the incidence of POAG was achieved by establishment of the competing endogenous RNA (ceRNA) regulatory network.¹¹ Through the integration of bioinformatics analysis and in vitro experimental investigations, the present study has revealed that ZFP42 exerts a mitigating effect on RGC damage in glaucoma by upregulating the expression of MARK2.¹² 

Mendelian randomization (MR) analysis, a statistical technique utilized in epidemiology to establish causal relationships, utilizes genetic variation as an instrumental variable.¹³ The utilization of MR analysis offers a reduced susceptibility to confounding, reverse causality, and measurement error as compared to conventional observational studies.¹⁴ MR analysis must satisfy the following three assumptions: (1) genetic variants selected as instrumental variables (IVs) are strongly relevant to exposure factors; (2) single nucleotide polymorphisms (SNPs) in genetic variants are not influenced by confounding factors that are associated with both exposure and result; and (3) genetic variants impact outcomes only through exposure and not through other biological mechanisms. IOP is a significant risk factor for the advancement of glaucoma; however, it is only detected in 70% of cases. Some patients, even with controlled IOP within the normal range, still show disease advancement.¹⁵ It is essential to investigate the risk factors for POAG to provide tailored screening for high-risk populations and prevent or postpone blindness in individuals. 

In this study, bioinformatics methods were employed to acquire expression matrices of optic nerve head (ONH) tissues from individuals diagnosed with POAG through the Gene Expression Omnibus (GEO) database, and differential expression analysis and enrichment analysis were conducted. Subsequently, MR analysis was integrated to investigate risk variables for POAG. Additionally, two machine-learning algorithms are discussed in this paper to identify biomarkers associated with POAG, and the potential mechanisms underlying their involvement in POAG were investigated by combining enrichment analysis, protein–protein interaction (PPI) networks, and subcellular localization analysis. This study offers novel perspectives on the genetic pathways that contribute to the initiation and advancement of POAG, along with possible biomarkers that could be utilized for tailored therapeutic interventions. 

Materials and Methods

Microarray Data Acquisition

The GSE45570—[HuEx-1_0-st] Affymetrix Human Exon 1.0 ST Array [transcript (gene) version]—was acquired from the GEO database (https://ww.ncbinlm.nih.gov/), containing six ONH samples from patients with POAG and six normal ONH samples. 

Identifying of Differentially Expressed Genes Between POAG and Control Groups

The differentially expressed genes (DEGs) between the POAG and control groups were identified via the limma R package (version 3.54.1; R Foundation for Statistical Computing, Vienna, Austria) with |log₂ fold change (FC)| > 0 and P < 0.05.^16,17 Additionally, the STRING database (http://string.embl.de/) was utilized to create a PPI network for these proteins so we could investigate protein–protein interactions. 

Enrichment Analysis

To comprehend the biological activities of these DEGs, we conducted an enrichment analysis via the clusterProfiler R package (version 4.6.0) (P < 0.05), based on the Gene Ontology (GO) database, including biological processes (BPs), molecular functions (MFs), and cellular components (CCs).¹⁸ 

MR Analysis

Dates of light stimulation, methylation, and POAG were gained from the IEU OpenGWAS database (https://gwas.mrcieu.ac.uk/). The childhood sunburn occasion (ukb-b-13246; 346,955 samples and 9,851,867 SNPs), frequency of solarium or sunlamp use (ukb-b-4943; 436,185 samples and 9,851,867 SNPs), sunburns (ebi-a-GCST90029034; 350,232 samples and 11,976,212 SNPs), and use of sun/ultraviolet protection (ukb-b-7422; 459,416 samples and 9,851,867 SNPs) were correlated with light stimulation. Moreover, six datasets were related to DNA methylation extraction, including Hannum-based DNA methylation age acceleration (ebi-a-GCST90014289; 34,449 samples and 7,541,726 SNPs); DNA methylation-estimated granulocyte proportions (ebi-a-GCST90014293; 6,152 samples and 12,859,948 SNPs); DNA methylation GrimAge acceleration (ebi-a-GCST90014294, 6148 samples and 12,962,189 SNPs; ebi-a-GCST90014300, 34,467 samples and 6,196,605 SNPs); Hannum-based DNA methylation age acceleration (ebi-a-GCST90014301; 6148 samples and 6,196,386 SNPs); and DNA methylation-estimated plasminogen activator inhibitor-1 levels (ebi-a-GCST90014303; 6149 samples and 6,196,880 SNPs). Also, the POAG-related data (bbj-a-75) contained 22,795 samples and 5961,428 SNPs. To satisfy the three basic assumptions of MR analysis, the SNPs used for analysis have to satisfy the following conditions: (1) SNPs are correlated with exposing factors (light stimulation or DNA methylation) at the genome-wide significance level (P < 5×10^–8); and (2) SNPs in linkage disequilibrium (r² = 0.001 and kb = 10,000) are excluded to ensure independence. 

MR analysis was carried out via five methods: weighted median, MR–Egger, simple mode, inverse variance weighted (IVW), and weighted mode. The results referred to are mainly IVW, as the IVW method provides more precise causation and its results are unbiased. Moreover, odds ratios (ORs) were calculated, with an OR value greater than 1 indicating that an exposure factor was a risk factor for the outcome, and an OR value less than 1 suggesting that an exposure factor inhibited the occurrence of the outcome event. The results are presented by scatterplots, forest plots, and funnel plots. 

To enhance the assessment of the accuracy of the findings, a sensitivity analysis was performed, which involved doing tests for heterogeneity and horizontal pleiotropy and employing the leave-one-out (LOO) method. In the heterogeneity test, a Q value greater than 0.05 indicates that there is no heterogeneity. P was greater than 0.05 in the horizontal pleiotropy test, suggesting that there was no horizontal pleiotropy. The LOO method was employed to detect any outliers with regard to the impact of each SNP. 

Screening of Biomarkers Related to Light Stimulation for POAG

Least absolute shrinkage and selection operator (LASSO; glmnet R package, version 4.1.4) and support vector machine-recursive feature elimination (SVM-RFE; e1071 R package, version 1.7-12) were utilized separately to identify the characteristic genes.^19,20 After that, the characteristic genes obtained by two machine-learning algorithms were intersected to obtain intersecting genes, which were noted as biomarkers. Furthermore, the diagnostic values of biomarkers for POAG were measured by drawing receiver operating characteristic (ROC) curves using the pROC R package (version 1.18.0).²¹ In addition, the expression levels of the biomarkers were evaluated between the POAG and control groups via the Wilcox test (P < 0.05). 

Subcellular Localization Analysis

To explore the specific location of biomarkers within cells, subcellular localization analysis was performed via the GeneCards database (https://www.genecards.org/). 

Gene Set Enrichment Analysis

We carried out enrichment analyses to further investigate which biological functions or signaling pathways are connected to biomarkers. The correlations between biomarkers and other genes were calculated, and all genes were ranked according to their correlation coefficients. Gene set enrichment analysis (GSEA) was performed via the clusterProfiler R package (version 4.6.0). 

Creation of the ceRNA Regulator Network

The miRNet database (https://www.mirnet.ca/) was utilized to predict the miRNAs regulated biomarkers. After that, the lncRNAs with regulatory interactions with miRNAs were predicted via the miRTarBase database (https://mirtarbase.cuhk.edu.cn/) and NPInter database (http://bigdata.ibp.ac.cn/npinter4/), and lncRNAs predicted by two databases were intersected to obtain common lncRNAs to conduct the ceRNA regulator network. Cytoscape software was utilized to visualize the ceRNA regulator network. 

Prediction of Small-Molecule Compounds

We employed the comparative toxicogenomics database (CTD; https://ctdbase.org/) to predict the small-molecule compounds correlated with biomarkers. The Cytoscape software was utilized to visualize the small-molecule compound–mRNA network. 

Correlation Between Biomarkers and DNA Methylation

First, DNA methylation–related genes were extracted from the GeneCards database (https://www.genecards.org/), including 64 double-strand break (DSB)-related genes, 106 CpG-related genes, and 53 m5c-related genes. Subsequently, these three types of genes were intersected to obtain common genes, and the Spearman algorithm was utilized to measure the relevance between biomarkers and these genes. Next, prediction of island CpGs for each biomarker was performed in a promoter region 2000 bp upstream via Methprimer (http://www.urogene.org/cgi-bin/methprimer/methprimer.cgi) and EBI CpGplot (http://www.ebi.ac.uk/Tools/seqstats/emboss_cpgplot/). Finally, a PPI network was established to investigate the interaction between light stimulation–related genes and DNA methylation–related genes. 

Construction of the Animal Model and Experiment Design

BALB/cJ mice, 2 to 5 months old, were obtained from the Beijing Charles River Laboratories, and mice were housed in a controlled environment with regulated temperature and humidity and a 12-hour light/dark cycle (lights on at 6:00 AM). Food and water were freely accessible. All animals were treated in line with the ARVO Statement for the Use of Animals in Ophthalmic and Vision Research, and the experimental protocol was approved by the Animal Ethics and Management Committee of Bengbu Medical College (approval no. 2023452). 

It has been shown that intravitreal injection of transforming growth factor-beta 2 (TGF-β2) in mice induces an increase in IOP.²² Therefore, in our study, all mice were stochastically divided into three groups, where one group of mice was injected in the eye with an adenoviral vector encoding a bioactive mutant of human TGF-β2 (Ad.hTGF-β2^226/228) to induce an increase in IOP. The other two groups included a control group and a group where the mice were injected in the eye with a null adenoviral vector (Ad.Empty). Before the vector was injected, mice were anesthetized with an intraperitoneal injection of a solution containing acepromazine (1.8 mg/kg), ketamine (73 mg/kg), and xylazine (1.8 mg/kg). One eye was randomly designated as the surgical eye for injection with one or two drops of 1% cyclopentolate (Mydriacyl; Alcon Laboratories, Fort Worth, TX) to dilate the pupil and one or two drops of 0.5% proparacaine (Alcaine; Alcon Laboratories) to serve as a local anesthetic. As previously described, 2 µL Ad.hTGF-β2^226/228 (6 × 107 pfu) or Ad.Empty (6 × 107 pfu) was injected intravitreally.²² The IOP was measured from 2 PM to 4 PM every day for 30 consecutive days. 

Reverse Transcription–Quantitative Polymerase Chain Reaction

The total RNA of mouse optic nerves was extracted using TRIzol (Invitrogen, Waltham, MA) in accordance with the instructions provided by the manufacturer. The QuantiTect Reverse Transcription Kit (Qiagen, Hilden, Germany) was utilized to perform reverse transcription of total RNA to cDNA. Reverse transcription–quantitative polymerase chain reaction (RT-qPCR) was carried out utilizing the SYBR Green Master Mix (Qiagen). The internal reference gene was β-actin. The primer sequences were β-actin (forward: 5′-CCATGTACGTAGCCATCC-3′; reverse: 5′-TCAGCTGTGGTGGTGAA-3′) and RAB8A (forward: 5′-GGATTCGGAACATTGAAGAGCA-3′; reverse: 5′-TCGAGTGCCAGCTTTTCTCC-3′). 

Statistical Analysis

R and Prism (GraphPad, Boston, MA) were utilized to perform the statistical analyses. The Wilcox test was employed to assess the differences between groups. P < 0.05 indicates a statistically significant difference. 

Results

Enrichment Pathways of DEGs are Associated With Light Stimulation and DNA Methylation

In total, 306 DEGs between the POAG and control groups were identified, containing 185 upregulated genes and 121 downregulated genes (Figs. 1A, 1B). The PPI network of DEGs that was created included 157 nodes and 232 edges (Fig. 1C). Angiotensinogen (AGT) had the strongest interactions with other proteins, as well as interactions with 13 proteins, such as melanin concentrating hormone receptor 1 (MCHR1), human glutamyl aminopeptidase (ENPEP), and neurotrophic receptor tyrosine kinase 2 (NTRK2). 

Figure 1.

 

Identification and enrichment analysis of DEGs comparing the POAG and control groups. (A) Volcano plot of DEGs. (B) Heatmap of DEGs. (C) PPI network maps of DEGs. (D) GO analysis of the top 10 functionally enriched DEGs. The red dots indicate upregulated genes, and the blue dots indicate downregulated genes. (E) The DEGs were related to light stimulation and methylation.

GO analysis revealed that 65 GO items were enriched, including 43 GO BP items, 12 GO CC items, and 10 GO MF items (Supplementary Table S1). The top 10 items of each type were identified (Fig. 1D). For example, of the GO BP items, the DEGs were enriched in visual perception, sensory perception of light stimulus, DNA demethylation, and so on. Of the GO CC, the DEGs were involved in photoreceptor outer segment, photoreceptor cell cilium, non-motile cilium, and so on. For GO MF, extracellular matrix structural constituent conferring compression resistance, photoreceptor activity, G-protein–coupled photoreceptor activity, and seven other items were enriched. In general, the majority of items that DEGs enriched were relevant to light stimulation and DNA methylation. Of all GO items, we screened the 25 GO items that were associated with light stimulation and DNA methylation (Supplementary Table S2). After de-embedding and merging, a total of 31 light stimulation–related genes (ASPN, B9D1, CRB1, CRYGS, GRK7, HTR2A, IMPG1, IMPG2, GUCA1B, MCHR1, MFAP5, NDP, NTRK2, OPN1SW, CDS1, EYS, RAB8A, RBP3, GNAT2, PRELP, SLC4A5, CALB1, OPN1LW, PRG3, GUCA1C, ARR3, CNGB3, SMAD3, PDE6H, NEUROD1, and ISL1) and five DNA demethylation–related genes (TDG, DPPA3, APOBEC3B, APOBEC3D, and APOBEC3F) were obtained (Fig. 1E). 

Three Light Stimulation–Related Biomarkers Were Obtained

Based on light stimulation–related genes, four characteristic genes were identified by the LASSO algorithm: RAB8A, PRG3, SMAD3, and OPN1LW (Figs. 2A, 2B). SVM-RFE identified a total of five characteristic genes: RAB8A, CNGB3, PRG3, SMAD3, and PRELP (Figs. 2C, 2D). Through intersecting, three biomarkers were obtained: RAB8A, PRG3, and SMAD3 (Fig. 2E); moreover, these three genes were all significantly highly expressed in the POAG group (Fig. 2F). Additionally, the area under the ROC curve (AUC) values of the three genes were 0.9722 (RAB8A), 0.8889 (PRG3), and 0.8889 (SMAD3), demonstrating that these genes had better diagnostic values for POAG (Fig. 2G). 

Figure 2.

 

Identification of biomarkers related to light stimulation. (A) LASSO regression coefficient path diagram. (B) LASSO cross-validation curves. (C) The curve of change for the predicted true values in the SVM-RFE algorithm. (D) The curve of change for the predicted error values in the SVM-RFE algorithm. (E) Venn diagram demonstrates the intersection of biomarkers obtained from the two algorithms. (F) The expression of biomarkers in patients with POAG compared with the control group. (G) ROC analysis of biomarkers.

Subcellular Localization of Biomarkers

Based on subcellular localization analysis, PRG3 was mainly extracellular (Fig. 3A). RAB8A was mainly located in the Golgi apparatus, cytosol, endosome, and cytoskeleton (Fig. 3B). Additionally, SMAD3 was mainly located in the cytosol and nucleus (Fig. 3C). 

Figure 3.

 

Subcellular localization analysis of biomarkers: (A) PRG3, (B) RAB8A, and (C) SMAD3.

Functional Enrichment Analysis of Biomarkers

The GO results indicated the top 10 items that were enriched. PRG3 was enriched in Golgi vesicle transport, transport along microtubule, axo-dendritic transport, etc. (Fig. 4A). RAB8A and SMAD3 were both involved in photoreceptor cell cilium, photoreceptor inner segment, detection of light stimulus, photoreceptor outer segment, visual perception, 9 + 0 non-motile cilium, and non-motile cilium, among others (Figs. 4B, 4C). 

Figure 4.

 

The GSEA of biomarkers. The GSEA in the GO analysis of PRG3 (A), RAB8A (B), and SMAD3 (C). The GSEA in the KEGG analysis of PRG3 (D), RAB8A (E), and SMAD3 (F).

The Kyoto Encyclopedia of Genes and Genomes (KEGG) results also identified enriched pathways. Protein processing in endoplasmic reticulum, olfactory transduction, ubiquitin-mediated proteolysis, and 73 other signaling pathways were related to PRG3 (Fig. 4D). RAB8A and SMAD3 were significantly negatively relevant to phototransduction and were markedly positively associated with human T-cell leukemia virus 1 infection and Epstein–Barr virus infection (Figs. 4E, 4F). In addition, RAB8A was markedly related to antigen processing and presentation, synaptic vesicle cycle, and influenza A, among others (Fig. 4E). SMAD3 was significantly correlated with the nuclear factor kappa B (NF-κB) signaling pathway, Kaposi sarcoma-associated herpesvirus infection, tumor necrosis factor (TNF) signaling pathway, etc. (Fig. 4F). 

Creation of a ceRNA Network and Small Molecule Compound–mRNAs network

Through prediction, eight miRNAs and 77 lncRNAs were obtained. The ceRNA regulator network contained three mRNAs, eight miRNAs, 77 lncRNAs, and 93 edges (including 16 miRNA–mRNA pairs and 77 miRNA–lncRNA pairs) (Fig. 5A). For example, the regulatory relationships included PRKAB1-hsa-mir-1224-3p-RAB8A, UHRF1-hsa-mir-124-3p-SMAD3, and DNAJA1-hsa-mir-335-5p-PRG3. Moreover, a total of seven miRNAs could regulate both RAB8A and SMAD3—namely, hsa-mir-1224-3p, hsa-mir-124-3p, hsa-let-7b-5p, hsa-mir-16-5p, hsa-mir-1-3p, hsa-mir-155-5p, and hsa-mir-23b-3p. In addition, hsa-mir-335-5p could regulate SMAD3 and PRG3. 

Figure 5.

 

Regulatory networks of biomarkers. (A) The ceRNA regulatory network; blue squares represent miRNAs, and yellow triangles represent lncRNAs. (B) A small molecule compound–mRNA network.

Based on the CTD, a small-molecule compound–biomarker network was created containing 66 nodes and 76 edges (Fig. 5B). Of these small-molecule compounds, benzo(a)pyrene and benzo(e)pyrene were associated with all of the biomarkers. 

DNA Methylation GrimAge Acceleration was Risk Factor for POAG

The IVW method indicated that DNA methylation GrimAge acceleration had a causal relationship to POAG (ebi-a-GCST90014294: P = 0.018, OR = 1.056; ebi-a-GCST90014300: P = 0.025, OR = 1.256) (Table) and was a risk factor for POAG. The results of the scatterplots showed that the slopes of the lines of DNA methylation GrimAge acceleration were positive, further validating the results (Supplementary Figs. S1A, S1B). Among the forest plot results, the overall effect values of the exposure factor on the outcome were greater than 0, which further supported the results (Supplementary Figs. S2A, S2B). The funnel plots demonstrated that MR adhered to Mendel's second law of random grouping (Supplementary Figs. S3A, S3B). 

Table.

 

IVW Results for Light Stimulation–Related Exposure Factors and Methylation-Related Exposure Factors

 

For the evaluation of the reliability of the results, we then performed sensitivity analyses. First, the Q values and P values were both greater than 0.05, suggesting that there was no heterogeneity and no confounding factors in this study (Supplementary Tables S3, S4). Thereafter, the LOO method suggested that there were no points of deviation (Supplementary Figs. S4A, S4B). In conclusion, DNA methylation GrimAge acceleration was a risk factor for POAG, suggesting that DNA methylation was essential in the development of POAG. 

Methylation in RAB8A Promoter Plays an Essential Role in POAG

The results for the prediction of CpG islands suggest that two CpG islands were found, at 400 to 700 bp and at 900 to 2000 bp, around the RAB8A promoter (Fig. 6A). For SMAD3, a CpG island was found at 1800 to 2000 bp around the promoter (Fig. 6B). However, for PRG3, no CpG island was predicted in the promoter (Fig. 6C). Therefore, these results demonstrate that DNA methylation of the RAB8A promoter might play an essential role in POAG. 

Figure 6.

 

Methylation analysis. (A–C) Prediction of CpG islands in RAB8A promoter, SMAD3 promoter, and PRG3 promoter. (D) Venn diagram demonstrates the intersection of DSB-related genes, CpG-related genes, and m5c-related genes. (E) The correlational heatmaps between biomarkers and DNA methylation-related genes. (F) Correlations among RAB8A and APOBEC3C, APOBEC3G, and APOBEC3D. (G) PPI network maps among light stimulation-related genes and APOBEC3C, APOBEC3G, and APOBEC3D. *P < 0.05, **P < 0.01, ***P < 0.001.

A total of 34 DNA methylation-related genes were then obtained, and correlation analysis showed that RAB8A was considerably associated with APOBEC3C (P = 0.01, r = 0.72), APOBEC3G (P = 0.04, r = 0.58), and APOBEC3D (P = 0.01, r = 0.73) (Figs. 6D–6F).The results of the PPI network further showed that low expression of HSP900AB1 could promote the expression of RAB8A and DNA methylation-related genes (including APOBEC3C, APOBEC3G, and APOBEC3D) (Fig. 6G). 

IOP was Significantly Induced and RAB8A was Markedly Highly Expressed in the Ad.hTGF-β2^226/228 Group

The results suggest that there was no significant difference in IOP between the control group and the Ad.Empty group (Fig. 7A). However, the IOP of mice in the Ad.hTGF-β2^226/228 group was markedly higher than that of mice in the other two groups (Fig. 7A). RT-qPCR demonstrated that the expression of RAB8A was highest in the Ad.hTGF-β2^226/228 group, and the difference was marked (Fig. 7B). This result was consistent with the above bioinformatics results. 

Figure 7.

 

Construction of animal models. (A) Changes of IOP with time in the control group, Ad.hTGF-β2226/228 group, and Ad.Empty group. (B) Expression of RAB8A in the control group, Ad.hTGF-β2226/228 group, and Ad.Empty group. *P < 0.05, **P < 0.01, ***P < 0.001.

Discussion

In this study, three biomarkers (SMAD3, PRG3, and RAB8A) associated with light stimulation were identified by combining bioinformatics and machine learning. Thereafter, the MR analysis demonstrated a causal relationship between DNA methylation (DNA methylation GrimAge acceleration) and POAG. Additionally, the correlation analysis further revealed that methylation of the RAB8A promoter may play a crucial role in the occurrence and development of POAG. 

Increased IOP significantly increases the likelihood of developing and worsening POAG, and heightened production of aqueous humor or blockage in its drainage can result in increased IOP.²³ Mouse models have established that the expression level of TGF-β2 is increased in the aqueous humor of POAG patients, and Smad3 is a crucial signaling protein for the TGF-β2–induced elevation of IOP and deposition of fibronectin in the trabecular meshwork.^24,25 In our study, SMAD3 expression levels were elevated in POAG, which is consistent with previous studies. Additionally, enrichment analysis revealed that SMAD3 had associations with the NF-κB signaling pathway, TNF signaling pathway, and other pathways. Saccà et al.²⁶ proposed that oxidation causes damage to the trabecular meshwork, resulting in major protein alterations in the aqueous humor. These changes may then lead to the diffusion of proteins from the front to the back of the eye, potentially triggering apoptosis of the RGCs.²⁶ The production of several pro-inflammatory cytokines, including TNF-α, interleukin 1 (IL-1), and IL-6, is triggered by reactive oxygen species–induced NF-κB activation, which intensifies the inflammatory response.²⁷ Numerous investigations have demonstrated that neuroinflammation plays a crucial role in the loss of RGCs in glaucoma, and TNF-α is believed to be a key contributor to this process.^28,29 In the aqueous humor of patients with glaucoma, the expression levels of TNF cytokines are higher than healthy levels.³⁰ Li et al.³¹ discovered that the activation of TRPV4 triggers gliosis in Müller cells and induces an increase in TNF-α levels through the JAK2/STAT3/NF-κB pathway, which has the potential to worsen RGC death in individuals with glaucoma. Our study demonstrated a positive relevance between the expression of SMAD3 and the NF-κB signaling pathway, as well as the TNF signaling pathway. Based on this, we hypothesized that SMAD3 would enhance the apoptosis of RGCs by activating the NF-κB signaling pathway and the TNF signaling pathway, consequently leading to the promotion of POAG. 

Glaucoma is classified as a neuroinflammatory condition due to the degeneration of RGGs and the infiltration of both active local glial cells and circulating immune cells.^32,33 In pathological conditions, neuroglia undergo a transformation into innate immune cells that function as antigen-presenting cells and release various chemicals, such as neurotrophic factors, into the retina to facilitate the healing of damage and eliminate severely damaged neurons.^34–37 Previous investigations have revealed that RRG3 increases neuronal protrusion and helps regulate filopodia development in immature neurons.³⁸ Moreover, RRG3 and PRG5 C-termini were found to play a crucial role in determining early neuronal morphology.³⁹ Su et al.⁴⁰ discovered that ARID1A influences the polarization of microglia during early fetal brain development by modifying the chromatin structure of microglial cells. Disturbances in microglial balance affect the release of PRG3, and the absence of PRG3 leads to changes in the downstream pathway of Wnt/β-catenin signaling by interacting with LDL receptor-related protein 6 (LRP6), a receptor on neural progenitor cells. These alterations result in disrupted regulation of neuronal development and the manifestation of autism-like behaviors in later stages of brain development.⁴⁰ Nevertheless, the role of PRG3 in glaucoma remains unexplored, and additional comprehensive investigations are required to ascertain the potential relevance between PRG3 and glial cells in the advancement of POAG. 

Previous research has demonstrated a correlation between mutations in proteins found in the optic nerve and the onset of glaucoma. Genetic evidence suggests that E50K (Glu50Lys) is a dominantly pathogenic mutation in photoreceptor protein that induces defective recycling of transferrin receptors in vivo through Rab8 inactivation mediated by the GTPase-activating protein TBC1D17.⁴¹ However, as of now, there are limitations in the research on the role of RAB8A in the development of glaucoma. 

Additionally, the MR analysis confirmed a causal connection between DNA methylation (DNA methylation GrimAge acceleration) and POAG. Previous research has shown that DNA methylation is involved in the progression of glaucoma. Human primary glaucomatous trabecular meshwork cells have elevated levels of DNA methylation, heightened expression of TGF-β1, and reduced expression of RAS protein activator like 1 (RASAL1) compared to normal trabecular meshwork cells.^42,43 Genome-wide DNA methylation studies of human primary Schlemm's canal endothelial cells demonstrated that POAG had different methylated regions (DMRs) of glaucoma-related genes when compared to the normals, further suggesting that DNA methylation is involved in the development of glaucoma.^44,45 In a study conducted by Wan et al., genome-wide DNA methylation microarray analysis of POAG and normal human trabecular meshwork samples showed that localized accumulation of GDF7 protein (a member of the TGF-β superfamily) promotes trabecular meshwork fibrosis, which in turn causes reduced aqueous outflow and increased IOP.⁴⁶ Furthermore, there is a significant association between DNA methylation and age, as evidenced by the observation that DNA methylation in juvenile retinas increases following several instances of elevated IOP.⁴⁷ These studies provide additional validation for the dependability of our study. In addition, the GpC island region was discovered to be highly prevalent in the RAB8A promoter in our study, suggesting that RAB8A is more prone to DNA methylation. The results of enrichment analysis suggest that RAB8A is relevant to the TGF-β signaling pathway. Prior research has indicated a potential connection among DNA methylation, TGF-β1, and RASAL1, which may contribute to the advancement of profibrotic diseases in the trabecular meshwork of glaucoma.⁴³ The findings from the PPI network suggest that the HSP90AB1–APOBEC3–RAB8A axis plays a significant role in POAG. Consequently, we propose that RAB8A might be linked to DNA methylation and the TGF-β signaling pathway in the development and advancement of POAG. 

The present investigation is subject to several constraints. The dataset GSE45570 used in this study has a small sample size, and further sequencing analysis with a larger sample size is required to validate the findings. Furthermore, this study lacked functional cellular and animal investigations to investigate the precise processes by which biomarkers contribute to the development of glaucoma. 

In summary, our research indicates that SMAD3, PRG3, and RAB8A are biomarkers for POAG that may play a role in the progression of the disease through neuroinflammation and fibrosis. Our research advances the study of POAG mechanisms and provides new targets for the treatment of POAG in the clinic. 

Acknowledgments

Supported by grants from the Anhui Engineering Technology Research Center of Biochemical Pharmaceutical (Bengbu Medical University, 22SYKFD01); Foundation of Educational Department in Anhui Province (KJ2021A0766, 2022AH051433); Foundation of Bengbu Medical University (2021byzd051, 2022byfy005, 2022byzd062); Foundation of Graduate Research and Innovation of Bengbu Medical College (Byycx23118); and China Prlmary Health Care Foundation. 

Disclosure: S. Zhao, None; Q. Dai, None; Z. Rao, None; J. Li, None; A. Wang, None; Z. Gao, None; Y. Fan, None