Potential Functions and Causal Associations of GNLY in Primary Open-Angle Glaucoma: Integration of Blood-Derived Proteome, Transcriptome, and Experimental Verification

Introduction

Primary open-angle glaucoma (POAG) is the predominant form of glaucoma characterized by progressive optic neuropathy, resulting in irreversible vision loss.¹ While increased intraocular pressure (IOP) remains the main risk factor, the pathogenesis of POAG is multifactorial, involving genetic predispositions, immune dysregulation, oxidative stress, and vascular dysfunction. Recent advances in genomic studies have significantly expanded our understanding of the genetic basis of POAG, particularly through genome-wide association studies (GWAS), which have identified hundreds of risk loci.^2,3 Despite these advances, translating genetic discoveries into mechanistic insights and clinical interventions remains challenging, as the causal genes and pathways underlying POAG susceptibility are often unclear.

The identification of genetic variants has greatly enhanced our understanding of glaucoma, a complex and heterogeneous group of eye disorders that cause optic nerve damage and vision loss. Research has pinpointed key genetic contributors to early-onset glaucoma, including the MYOC, OPTN, and TBK1 genes, which play pivotal roles in genetic testing and counseling.⁴ For primary open-angle glaucoma (POAG), over 74 genomic regions have been linked to disease susceptibility, highlighting the genetic complexity and variability underlying glaucoma pathogenesis.⁵ Specific subtypes, such as normal-tension glaucoma (NTG), characterized by optic nerve damage without elevated intraocular pressure, are associated with mutations in genes like OPTN and TBK1, with genetic predisposition observed in up to 21% of cases.⁶ Advances in genome-wide association studies (GWAS) have further expanded our understanding of glaucoma genetics.^2,3 These studies have identified more than 20 genetic loci associated with various glaucoma forms, illuminating novel biological pathways involved in disease progression.⁷ GWAS findings also reveal genetic overlaps between glaucoma and other ocular traits, including intraocular pressure regulation and vertical cup-to-disc ratio variations.⁵ Polygenic risk scores (PRS), derived from GWAS data, have emerged as powerful tools for assessing the genetic burden of POAG. These scores facilitate predictions regarding disease onset, intraocular pressure trends, and the likelihood of requiring surgical intervention. However, their applicability across diverse ethnic groups remains limited due to the predominant focus on Northern European populations.^8,9 Genetic testing for early-onset glaucoma genes provides substantial benefits for pre-symptomatic detection and patient counseling.⁴ Looking ahead, gene-based therapies hold great promise for personalized treatment strategies, targeting specific genetic pathways to meet individual patient needs. Such advancements aim to address current gaps in therapeutic options, enabling more precise and effective glaucoma management.^5,10

Extensive data on quantitative trait loci (QTL) have been published recently to show connections between gene expression (eQTL) and genotype and protein abundance (pQTL).^11,12 This progress has spurred the development of various statistical methods focused on efficiently integrating multidimensional data.¹³ Genes associated with protein levels in depression and Alzheimer’s disease have been found using proteome-wide association studies (PWASs).^14,15 Likewise, transcriptome-wide association studies (TWASs) have been conducted to examine the connections between different phenotypes and gene expression.¹⁶ In addition, Mendelian randomization (MR) techniques have been extensively used to pinpoint potential genes via the combination of QTL and disease-related GWAS data.^17,18 Based on three fundamental presumptions, MR may provide light on causation by simulating a natural randomized controlled trial in which mutations are allocated at random during meiosis.^19,20 In general, the amalgamation of GWAS data with multidimensional QTL data has the capacity for pinpointing certain pathways and potential genes that are crucial for understanding the fundamental workings of glaucoma.

Moreover, an earlier GWAS study using Mendelian randomization and genetic correlation-based approaches has also reported an association between glaucoma and immune-related diseases.² Studies have shown that T-cell responses, particularly those targeting heat shock proteins (HSPs), play a significant role in the degeneration of nerves in glaucoma. Elevated levels of HSP-specific T helper type 1 (Th1) cells and corresponding cytokines, such as interferon-γ (IFN-γ), have been observed in patients with POAG, correlating with thinner retinal nerve fiber layer thickness (RNFLT).²¹ Both experimental models and human studies suggest glaucoma may be an autoimmune disease, with antibodies and T-cell responses to HSPs detectable in glaucomatous patients and animal models.²² Analysis of T-cell subsets in glaucoma patients reveals a shift towards a proinflammatory state, characterized by decreased regulatory T cells (Tregs) and increased proliferative activity and cytokine secretion in CD4⁺ T cells.²³ Additionally, Th17 cells, known for promoting autoimmunity, are implicated, with glaucoma patients exhibiting higher levels of Th17-promoting cytokines and an enhanced capacity of Th17 cells to stimulate immunoglobulin production.²⁴ We therefore used our investigation on blood to get a more thorough grasp of the pathophysiology of glaucoma.

In this study, we aim to bridge the gap between genetic discoveries and disease mechanisms by leveraging integrative multi-omics approaches. Utilizing GWAS datasets, coupled with pQTL and eQTL data, we conducted PWAS, TWAS, and summary data-based Mendelian randomization (SMR) analyses to identify genes with causal roles in POAG. Our comprehensive analysis pinpointed granulysin (GNLY) as a key gene, supported by both genomic data and experimental validation. GNLY and its interacting genes were further analyzed for their roles in immune modulation, apoptosis, and cytotoxicity, highlighting their potential contributions to the neurodegenerative processes in POAG. This work establishes a framework for dissecting the genetic underpinnings of POAG and paves the way for novel therapeutic strategies. See Figure 1 for a flow chart depicting the study selection process.

Figure 1 (A) Shows the study design as a text-based flowchart. (B) Presents the same study design as a graphical flowchart.

Materials and Methods

Glaucoma GWAS Data

The glaucoma GWAS summary datasets (https://gwas.mrcieu.ac.uk/datasets/ebi-a-GCST90018852/) used in this investigation were taken from the FinnGen R10.²⁵ The main study of the glaucoma GWAS summary datasets was limited to people of European heritage, including information on 484,979 individuals (10,411 primary open-angle glaucoma cases and 474,568 controls) (Table S1). This large sample size enables robust analysis and enhances the statistical power to detect significant associations. The GWAS data underwent rigorous quality control (QC) procedures, including assessments of genotyping accuracy, minor allele frequency, and adherence to Hardy–Weinberg equilibrium. Following these QC measures, the final dataset included a specific number of single nucleotide polymorphisms (SNPs) that met the established criteria. The exact number of SNPs will be reported upon completion of the QC process.²⁵

Data on Human Blood Proteins and Transcripts

A major demographically determined research (AGES Reykjavik study; n=5457) provided the serum proteome data (Table S2).²⁶ The AGES Reykjavik research included people over 65 who were mostly European and for whose phenotypic and genetic data were available. 4137 human proteins were quantified in serum utilising the proteome profiling platform Slow-Off rate Modified Aptamer (SOMAmer), and their weights were used for PWAS.²⁶

The Genotype-Tissue Expression (GTEx) version 8 database provided the whole-blood eQTL data (Table S2)²⁷ and their corresponding mRNA weights can be accessed at http://gusevlab.org/projects/fusion. Paired-end sequencing (Illumina TruSeq; Illumina Inc.) was used to gather the RNA-seq data, and whole-genome sequencing was employed to determine genotypes. Visit the official GTEx website for further details on donor registration, consent processes, sample collecting techniques, sample processing, and histopathological analysis.

GNLY-Interacting Genes Network and Gene Set Enrichment Analysis

The online Search Tool for the Retrieval of Interacting Genes²⁸ (STRING; 2024 release) database was used to build the network of GNLY-interacting genes in order to assess the interaction linkages between the genes. A cumulative score greater than 0.9 was used to classify interactions as statistically noteworthy. Version 3.9.1 of the Cytoscape programme²⁹ was used to visualise the included regulatory networks.

The “cluster Profiler” package in R software was used for Gene Ontology function (GO) and Kyoto Genome and Genome Encyclopedia (KEGG) enrichment analyses to elucidate the functions and pathways that were enriched in the GNLY-interacting genes associated with POAG, and factors with an adjusted p value < 0.05 were considered statistically significant. GO lollipop plots, and KEGG lollipop plots were drawn. The GO analysis focused on molecular functions, cellular components and biological processes, whereas the KEGG analyses identified enriched signaling pathways.

Collection of Human Blood Samples

This investigation included individuals aged 40 to 65 years with primary open-angle glaucoma (POAG), encompassing both newly diagnosed inpatients and outpatients. Clinical data were collected from a total of 16 POAG patients and 44 healthy controls. From this cohort, peripheral blood samples were randomly selected from 5 POAG patients and 5 healthy controls for enzyme-linked immunosorbent assay (ELISA) analysis of GNLY expression levels. This selection was made to balance experimental feasibility with cost considerations, while ensuring representation within the clinical dataset. Peripheral blood specimens were collected at the Ophthalmology Department of the First Affiliated Hospital of Chongqing Medical University. Exclusion criteria included significant ocular diseases, recent ocular surgery, autoimmune diseases, malignancy, diabetes mellitus, steroid use, and smoking, to eliminate potential confounding factors. All participants provided informed consent, and the study was conducted in compliance with ethical standards, with approval from the Ethics Committee at the hospital (2024–098-02). This methodological approach ensured a meaningful comparison of GNLY expression levels between POAG patients and healthy controls, while adhering to rigorous ethical and procedural standards. Serum GNLY concentrations were measured using Human Granulysin Elisa kit (Cat# RX105613H; Ruixinbio Quanzhou, China). The sera of POAG patients and controls were prepared from peripheral venous blood samples (3 mL/sample), which were taken by venipuncture in 5mL plastic serum tubes (KWS 240914), allowed to clot for 20 min at 22–25°C, and centrifuged (4000 rpm, 20 min, at 4°C). The supernatant (serum) was collected in 1.5 mL, round-bottom microcentrifuge tube (BS-15-MA; Biosharp) and stored at −80°C until required for ELISA. Assays were performed following the manufacturer’s instructions, and the absorbance was measured using an MRX Revelation microplate reader.

Cell Validation

Mouse retinal ganglion cells (RGCs) (iCell Bioscience Inc., Shanghai, China) were grown in a mixture of F12 medium (DMEM/F12, Gibco, USA) supplemented with 20% fetal bovine serum (FBS, Gibco, USA) and 1% penicillin/streptomycin (Procell, Wuhan, China), as described by Park et al.³⁰ They were identified by anti-rat Thy1.1 mAb (1:200, Abcam), and anti-rat Map2 pAb (1:200, Abcam). Specific fluorescent pictures were provided in the supplementary material (Supplementary Figure S1). The cells were cultivated in a regulated setting with a humidity level and temperature maintained at 37°C, and a CO2 concentration of 5%. Excitotoxicity was induced by modifying the culture conditions with 2.5 mM N-methyl-D-aspartic acid (NMDA) (MCE, Shanghai, China) for a duration of 8 hours, following the methodology outlined in earlier investigations.^31,32 The GNLY concentrations of mouse retinal ganglion cells were measured using the Mouse Granulysin Elisa kit (Cat# RX2D2014796; Ruixinbio Quanzhou, China). The mouse retinal ganglion cells of POAG and controls were were harvested by trypsinization and centrifuged (1000 x g, 5 min). The cells were collected in 1.5 mL, round-bottom microcentrifuge tube (BS-15-MA; Biosharp) and stored at −80°C until required for ELISA. Assays were performed following the manufacturer’s instructions, and the absorbance was measured using an MRX Revelation microplate reader.

Statistical Analysis

We employed the Functional Summary-based Imputation (FUSION) to conduct PWAS and TWAS, a robust approach that combines gene expression/protein abundance data with summary statistics from GWAS to detect genes exhibiting associations between their cis-regulated expression/protein abundance and complex traits. We implemented the FUSION pipeline using the default settings, following the investigator’s guidelines to conduct PWAS and TWAS. FUSION employs multiple predictive models (top1, susie, lasso, enet) to determine the collective impact of single nucleotide polymorphisms (SNPs) on the weights assigned to gene (http://gusevlab.org/projects/fusion/proteinexpression).²⁶ Furthermore, 1000 Genomes European samples (https://data.broadinstitute.org/alkesgroup/FUSION/LDREF.tar.bz2) were utilized to estimate the Linkage disequilibrium (LD) between the prediction model and the SNP at each locus of GWAS. Each gene underwent a permutation test, repeated 2000 times, where the significance of observed differences was assessed both through a Z-test (assuming normal distribution, p<0.05) and the permutation test (which does not rely on distribution assumptions, p<0.05), providing a robust analysis of gene expression differences. To comprehensively identify proteins/genes associated with primary open-angle glaucoma (the lowest p-value in the PWAS dataset was > 0.05), we considered a p-value threshold of <0.05 as statistically significant for subsequent analyses, a criterion that has been frequently employed in previous studies^33–35 (Tables S3 and S5).

The SMR software tool was developed to implement the SMR and HEIDI methods. It uses pooled data from GWAS and eQTL studies to analyze how gene expression levels are pleiotropically associated with complex traits due to shared and potentially causal variants at specific loci.³⁶ The SMR and HEIDI approaches analyze whether gene expression mediates the effects of SNPs on POAG. Therefore, this tool helps prioritize genes identified in GWAS for further functional studies. In SMR, the top cis-eQTL estimates the impact of gene expression on the outcome, and this analysis hinges on the reliability of genetic variants as instrumental variables for consistent estimates. A genetic variant must satisfy three assumptions to be considered a valid instrumental variable IV: 1) it must be associated with the risk factor (gene expression); 2) it must not be associated with any confounders of the risk factor-outcome relationship; and 3) it must be conditionally independent of the outcome, given the risk factor and confounders. It was found that a p-value of less than 1×10⁻⁵ in the linear regression of the eQTL analysis for each variant indicates a small weak instrument bias.³⁷ In our SMR analysis, we used a standard threshold of P_eQTL = 5×10⁻⁸ to identify the most significant associated cis-eQTL, with a default window size of 2000kb. This method aims to minimize the inclusion of weak instrumental variables. Additionally, we systematically excluded SNPs with an allele frequency difference greater than 0.2 in any pairwise dataset. This includes the linkage disequilibrium (LD) reference sample data, the eQTL summary data, and the GWAS summary data. We also conducted the heterogeneity in dependent instruments (HEIDI) test to evaluate the possible linkage within the observed association. Rejection of the null hypothesis (ie, P_HEIDI < 0.05) indicates that the observed association could be due to two distinct genetic variants in high linkage disequilibrium with each other. In conducting the HEIDI test, we adopted the default settings: removing SNPs in very strong linkage disequilibrium (LD, r² > 0.9) with the top associated eQTL, and SNPs in low LD or not in LD [r² < 0.05] with the top associated eQTL; P_eQTL < 3.7 × 10⁻⁸; number of cis-SNPs≥3 for a HEIDI as HEIDI test loses power if cis-SNPs< 3; and maximum eQTLs in a HEIDI test = 20. Using the command line interface (https://yanglab.westlake.edu.cn/software/smr/#Overview), SMR analysis was performed using the default parameters (Table S4).³⁶

We did not correct for multiple testing in this exploratory study. All statistical data analyses were two-sided at a 0.05 level of significance and conducted with R software. The SMR, and FUSION software were used to facilitate the SMR, TWAS, PWAS analyses and display the results.

Results

PWAS of Glaucoma

Within our research, we determined 86 genes (SFTPD, CSK, COL18A1, TCN2, GZMK, RAB2A, TEK, GNLY, etc.) via PWAS. The procedure consisted of combining protein data obtained from AGES with glaucoma GWAS results (Figure 2A and Table S5). The PWAS procedure utilizes the integration of genetic data and protein quantity data to discover genes that might impact the probability of glaucoma by controlling protein levels in the bloodstream. The importance of these genes was assessed using a strict p value threshold set at 0.05. Furthermore, we integrated the cis-eQTL summary data from GTEx V8 with glaucoma GWAS data for SMR, identifying 405 more genes and taking an intersection of SMR genes and PWAS genes to get eight genes (SFTPD, CSK, COL18A1, TCN2, GZMK, RAB2A, TEK, and GNLY) as candidate casual genes (Figure 2B and Table S6). We have discovered a total of eight potential risk genes for glaucoma by analyzing both pQTL and eQTL data.

Figure 2 The circular Manhattan plots display the glaucoma PWASs in the proteomes of human plasma. A Manhattan figure was created to visualize the integration of the glaucoma GWAS with the AGES proteins (n = 5457) (A). Each segment on the circular axis represents a chromosome, labeled 1 through 22 for autosomes. Each dot on every chromosomal segment signifies a gene, while the intensity of association on the radial location indicates the -log10(p) of PWAS. The significance threshold for the whole proteome in the AGES dataset was established at p < 0.05, which is where the black circle is located. Venn diagram for SMR genes and PWAS genes (B). Proteome-wide significant genes (SFTPD, CSK, COL18A1, TCN2, GZMK, RAB2A, TEK, and GNLY) both in results of the SMR analysis and PWAS analysis are shown in red in the circular Manhattan plots. Blue stands for SMR and red for PWAS.

TWAS of Glaucoma

Since proteins are synthesized from mRNA, we hypothesized that genes linked to primary open-angle glaucoma (POAG) at the protein level could also show comparable evidence at the transcriptional level. Therefore, TWAS was conducted using blood eQTL data. 429 genes (AMZ2, BCL7A, CACFD1, DPP3, FZD2, GNGT2, HPSE, KIF5B, etc)., were linked to POAG at the level of gene expression. This method included merging transcriptome data from GTEx V8 with glaucoma GWAS results (Figure 3A and Table S3). Moreover, we took an intersection of PWAS genes and TWAS genes to get 8 genes (SIGLEC12, FCGR2B, C4BPA, APOL1, HAVCR2, SMAP1, PTGDS, and GNLY) as candidate casual genes (Figure 3B and Table S7). In total, we identified 8 POAG-related genes (SIGLEC12, FCGR2B, C4BPA, APOL1, HAVCR2, SMAP1, PTGDS, and GNLY) detected by PWAS and TWAS analyses.

Figure 3 Manhattan plots for the glaucoma TWASs in the human whole blood transcriptome. Manhattan figure for the TWAS merging the glaucoma GWAS with the GTEx V8 transcriptome (n = 558) (A). The x-axis shows individual genes, with the degree of correlation shown by the -log10(p) of TWAS on the y-axis. In the GTEx V8 dataset, the threshold for transcriptome-wide significance was set at p < 0.05, which is where the red horizontal line is located. Venn diagram for PWAS genes and TWAS genes (B). Eight common proteins/genes (SIGLEC12, FCGR2B, C4BPA, APOL1, HAVCR2, SMAP1, PTGDS, and GNLY) both in results of the PWAS analysis and TWAS analysis are shown in black in the Manhattan plots. Blue stands for TWAS and red for PWAS.

SMR Verified Genes Having Causal Associations with Glaucoma

Utilizing the GTEx V8 cis-eQTL statistics, our SMR approach has effectively revealed 405 genes that exhibit a non-pleiotropic association with glaucoma. These genes, namely SEMA3A (ENSG00000075213; β[SE]= −1.50[0.37], P_SMR = 5.81×10⁻⁵), HIVEP3(ENSG00000127124; β[SE] = −0.54 [0.14], P_SMR = 1.62×10⁻⁴), and RAB2A (ENSG00000104388; β[SE] = −0.50 [0.14], P_SMR = 3.89×10⁻⁴), have been selected as the top three genes during our research (Table S4). The majority of the discovered genes are situated on chromosome 6, chromosome 2, chromosome 1, and so on. Additionally, we took an intersection of SMR genes, TWAS genes and PWAS genes to get one gene (GNLY; ENSG00000115523; β[SE] = −0.16 [0.08], P_SMR = 4.17×10⁻²) as a potential risk gene (Figure 4A and Table S8). The SMR analysis identifies a significant association between the genetic variants of the GNLY gene and glaucoma. The top cis-eQTL, depicted in the scatter plot, suggests a notable inverse relationship between eQTL and GWAS effect sizes, indicating that the genetic expression alterations in GNLY may influence the susceptibility to glaucoma (Figure 4B). Additionally, the Manhattan plot highlights a peak that surpasses the suggestive significance threshold (P_SMR = 0.05), reinforcing the possible function of GNLY genetic variations in the development of glaucoma (Figure 4C).

Figure 4 Venn diagram for SMR genes, TWAS genes and PWAS genes (A). Blue represents TWAS, red represents PWAS, and purple represents SMR. Association of eQTL and GWAS effect sizes for GNLY (B). It illustrates the relationship between the effect sizes of expression quantitative trait loci (eQTL) on the x-axis and genome-wide association study (GWAS) effect sizes on the y-axis for the gene GNLY (ENSG00000115523). Each point represents a specific eQTL, with blue circles denoting cis-eQTL, and the red triangle indicating the top cis-eQTL. Error bars along both axes represent confidence intervals for effect sizes. A dashed Orange line is included to suggest a potential linear relationship between the two variables. Pleiotropic relationship between POAG and GNLY (C). Top plot: Solid rhombuses indicate that the probes pass the HEIDI test, whereas grey dots show the −log10 (p values) for SNPs from the POAG GWAS. Results of eQTL, middle plot. Gene locations marked by the probes are shown in the bottom plot. SMR, summary data-based Mendelian randomization; HEIDI, heterogeneity in dependent instruments; GWAS, genome-wide association study; eQTL, expression quantitative trait loci; POAG, primary open-angle glaucoma.

GNLY-Interacting Genes Network and Gene Set Enrichment Analysis

Integrative analyses using PWAS, TWAS, and SMR identified GNLY as the key gene. Further exploration through the STRING database revealed ten additional genes interacting with GNLY, none of which were identified in the initial analyses. Subsequent GNLY-interacting genes network was visualized by Cytoscape software (Figure 5A). Functional enrichment analyses of these genes via KEGG and GO pathways provided insights into the biological processes and pathways potentially regulated by GNLY and its network partners.

Figure 5 Network diagram of GNLY-interacting genes. The circle size and color intensity represent the number of connections (degree) of each gene, the larger the circle and the darker the color indicate the more important it is in the network. The edges (lines) indicate the relationship between genes, and more connections indicate higher network correlation between genes (A). Lollipop plots generated by GO/KEGG enrichment analysis.BP, biological process; CC, cellular component; MF, molecular function. Lollipop length represents the adjusted P value in the enrichment pathway, and bubble size represents the number of genes (B).

The GO and KEGG enrichment analyses of the DEGs showed that the functions of the GNLY-interacting genes were mainly enriched in cell killing, leukocyte mediated cytotoxicity, immunological synapse, and serratase activity (Figure 5B). The Kyoto Genome and Genome Encyclopedia enrichment analysis showed that they were mainly enriched in Natural killer cell mediated cytotoxicity, and Allograft rejection signalling pathway (Figure 5B and Table S9).

Summary of Findings

Levels of confidence in potential POAG-related genes were assessed using results from PWAS, SMR, and TWAS. GNLY was identified as having a very high confidence level for its causal effect on POAG, with all three analyses supporting its causality. Additionally, high confidence was found for the causal effect of SFTPD, CSK, COL18A1, TCN2, GZMK, RAB2A, and TEK on POAG, as PWAS and SMR analyses supported their causality. However, low confidence was observed for SIGLEC12, FCGR2B, C4BPA, APOL1, HAVCR2, SMAP1, and PTGDS. While these genes were detected in PWAS and TWAS, SMR did not support their causality with POAG, suggesting they may be associated with the risk for POAG rather than being causal factors (Table 1).

Table 1 The Compilation of Genes That Have Been Shown to Be Consistent Across at Least Two Datasets/Methods

Validation of the Hub Gene GNLY by Clinical and Cellular Experiments

An ELISA assay was used to confirm the expression of GNLY in the POAG groups. The results showed a substantial decline in the protein levels of GNLY in the cell model categories in comparison to those in the control categories (Figure 6A). Consistent with the above, the level of GNLY expression in the circulating blood of individuals diagnosed with POAG was much lower compared to that seen in healthy individuals (Figure 6B). Analysis of clinical data revealed that patients with POAG exhibit compromised visual function, inflammatory processes and kidney function. These manifestations primarily manifest as visual field defects, reduced thickness of the retinal nerve fiber layer (RNFL), larger cup-to-disc ratio, increased intraocular pressure (IOP), increased lymphocyte count, monocyte count and neutrophil count, higher creatinine levels (Table 2 and Additional file 2).

Table 2 Comparative Analysis of Preliminary Clinical Data Comparing Patients with POAG and Those Without the Condition

Figure 6 Validation of protein expression of GNLY. (A) ELISA to detect protein expression levels of GNLY between POAG mouse cell model groups and control groups. (B) ELISA to detect protein expression levels of GNLY between POAG patients and control subjects. ***p<0.001, p value for Mann–Whitney U-test. ELISA, enzyme-linked immunosorbent assay; POAG, primary open-angle glaucoma. Black represents the control group and gray represents the POAG group.

Discussion

Glaucoma is an intricate and multifaceted neurological disease, where several factors such as biomechanics, blood flow, metabolism, oxidative stress, and inflammation all play significant roles in its pathophysiology.³⁸ In addition, epigenetic aging acceleration is strongly associated with the progression of glaucoma.³⁹ The accumulation of senescent cells in the trabecular meshwork and retinal ganglion cells (RGCs) contributes to glaucomatous damage.⁴⁰ Cellular senescence also impairs the retina’s ability to recover from stress, such as elevated intraocular pressure, reducing its adaptability and increasing susceptibility to glaucoma.⁴¹ Moreover, oxidative stress-induced shortening of retinal ganglion cell telomeres establishes a critical link between aging and the pathogenesis of glaucoma.⁴² Consistent with prior studies, our clinical observations revealed that individuals with POAG displayed more pronounced ocular anatomical and functional deficits in contrast with the controls. The impairments seen in the study included increased IOP, increased cup-to-disc ratio, decreased best-corrected visual acuity, visual field abnormalities, and reduced thickness of the retinal nerve fiber layer (RNFL). Our research disclosed a noticeable trend towards an increased lymphocyte count, monocyte count and neutrophil count in the group with POAG compared to the controls. Although glaucoma is common, there is a need to improve and develop current diagnostic and treatment methods. Using information from the blood, we conducted a comprehensive examination of PWAS, TWAS, and SMR in this inquiry.

We identified 86 candidate risk genes for glaucoma in the AGES datasets via PWAS. Each group of proteins found may be responding to distinct components of the disease process, indicating the intricate and multifaceted characteristics of glaucoma. According to SMR analysis, 8 of these genes (SFTPD, CSK, COL18A1, TCN2, GZMK, RAB2A, TEK, and GNLY) have a robust causal correlation with glaucoma. While, two of these (COL18A1, and TEK) have previously been described in correlation with glaucoma. The COL18A1 gene, encoding collagen type XVIII, is implicated in primary angle-closure glaucoma (PACG) through mutations that potentially cause iridocorneal angle closure, a critical feature of PACG.⁴³ Variations in COL18A1 are also associated with age of onset in primary open angle glaucoma (POAG).^44,45 Similarly, TEK is a receptor tyrosine kinase that plays a key role in the development of Schlemm’s canal and trabecular meshwork associated with primary congenital glaucoma (PCG).^46–48 These TEK mutations result in haploinsufficiency, reducing protein function and impacting cellular mechanisms like protein degradation and ligand responsiveness.^49,50 Additionally, the angiopoietin-TEK signaling pathway plays a significant role in vascular development and is implicated in the etiology of PCG.^48,51 However, other genes are newly identified, and their role in glaucoma has not been previously reported, warranting further inquiry.

Then, we identified eight genes, including SIGLEC12, FCGR2B, C4BPA, APOL1, HAVCR2, SMAP1, PTGDS, and GNLY, which were disclosed by both PWAS and TWAS. This demonstrates that the findings of this research exhibit coherence in both the process of converting genetic information into proteins (translation) and the process of copying DNA into RNA (transcription). However, excluding PTGDS, none of these have been previously documented in connection with glaucoma. Prostaglandin D2 synthase (PTGDS) expression has been linked to the regulation of intraocular pressure (IOP). Genetic variations related to prostaglandin pathways, specifically PTGS1 and PTGFR, have been associated with IOP response to latanoprost, a common glaucoma medication.⁵²

Interestingly, only GNLY was further supported by SMR. We have identified a substantial association between the expression of GNLY in whole blood and a lower risk of glaucoma. These findings indicate that the discovered genes may have an influence on the control of the development of glaucoma and might be suitable subjects for further investigation. The findings of our research emphasize the intricate nature of the genetic foundation of glaucoma. The GNLY gene encodes granulysin, a cytolytic antimicrobial peptide (AMP) crucial in the immune response against a variety of pathogens and in mediating inflammatory diseases. It has been shown that granulysin is significantly involved in the pathogenesis of psoriasis, where its expression is increased in lesions, and certain polymorphisms within the GNLY gene, such as rs7908, are associated with the disease’s severity.^53–55 Moreover, granulysin plays a pivotal role in the immune defense against viruses, including hepatitis B virus (HBV), suggesting that genetic variations in the GNLY gene could influence host resistance to HBV infection.⁵⁶ In addition to its role in psoriasis and viral infections, granulysin is also implicated in severe cutaneous adverse reactions like Stevens–Johnson syndrome (SJS) and toxic epidermal necrolysis (TEN). A study identified a GNLY genetic variant that contributes to the pathophysiology of SJS/TEN, highlighting the importance of granulysin in skin health and disease.⁵⁷ Furthermore, granulysin’s cytotoxic activity extends to its potential therapeutic application in cancer treatment, as demonstrated by enhanced tumor cell killing in GNLY transgenic mice.⁵⁸

The analysis of GNLY-interacting genes’ functional enrichment reveals their involvement in crucial biological processes like immune dysregulation, inflammation, and apoptosis. These processes are pivotal in the pathophysiology of primary open-angle glaucoma (POAG), particularly in the context of retinal ganglion cell (RGC) death, a hallmark of glaucoma. They are intricately connected with other pathogenic mechanisms such as elevated intraocular pressure (IOP), oxidative stress, and vascular dysfunction. Neuroinflammation plays a critical role in POAG, with immune and glial cell activation contributing to disease progression.^59–62 Abnormal immune responses, encompassing both protective and detrimental effects, are associated with optic nerve degeneration in glaucoma, emphasizing the importance of balanced immune regulation in mitigating neurodegeneration.^22,63–65 Moreover, oxidative stress, a significant trigger of immune dysregulation, is closely associated with RGC apoptosis. Pathways like NF-κB activation and pro-inflammatory cytokine signaling are commonly observed in glaucoma patients.⁶⁶ Apoptotic mechanisms, including Fas-mediated and Bcl-2-dependent pathways, are crucial drivers of RGC loss, with their regulation offering potential therapeutic opportunities.⁶⁷ Elevated IOP-induced mechanical stress further worsens neuronal apoptosis, interacting with factors such as neurotrophic factor deprivation and vascular dysregulation to contribute to disease progression.⁶⁸ Additionally, pathways involving NF-κB, complement activation, and toll-like receptors play critical roles in cellular injury and neuroprotection, emphasizing the multifaceted nature of POAG pathogenesis.⁶⁹ Collectively, these findings offer a comprehensive insight into the biological significance of GNLY-interacting genes in the context of POAG.

The diverse range of proteins used in this study are crucial in elucidating the potential biological processes associated with glaucoma. By conducting pQTL analysis, these proteins have been linked to genetic variability, investigating the influence of genetic mutations on blood protein levels. Nevertheless, our research is subject to many constraints. First, one of the limitations of this study is the relatively small sample size used for validating GNLY expression, which may limit the statistical power and generalizability of the findings. Future studies should include larger cohorts with diverse demographic and clinical characteristics to confirm these results and ensure their applicability across broader populations. Additionally, multicenter collaborations could help overcome sample size restrictions and provide a more comprehensive understanding of GNLY’s role in POAG. Second, the dataset size for glaucoma GWAS and pQTL information was restricted, potentially impacting the reliability of our results. With the increasing availability of glaucoma GWAS and pQTL datasets in the future, the importance and applicability of PWAS in comprehending disorders such as glaucoma are expected to become increasingly apparent. Then, the fact that our research was limited to European people may have affected the discovery of specific impacts of gene transcriptome and proteome expression. This restriction hinders the applicability of our results, underscoring the necessity for additional research involving more diverse populations to corroborate our findings. Finally, we only focused on a single tissue relevant to glaucoma, potentially limiting our exploration of connections between glaucoma and other tissues. It is possible that certain transcripts are found in the retina but not in the bloodstream. Of course, although the current study provides evidence of GNLY involvement in POAG through ELISA and in vitro experiments, these findings require in vivo validation to establish their physiological relevance. Future studies using suitable animal models, such as transgenic or knockout mice, could provide greater insights into the pathogenic mechanisms of GNLY and its interaction with related genes. Such models would also allow for the evaluation of potential therapeutic interventions targeting GNLY-mediated pathways in POAG. Notably, there is still insufficient evidence in new sequencing data, particularly at the single-cell and tissue levels. Additional investigation should concentrate on investigating the mechanisms through which GNLY influences the pathophysiology of glaucoma using experimental approaches. Although there were limitations, we successfully integrated these datasets to identify GNLY as a possible cause of glaucoma and acquire understanding of its possible mechanisms of action.

Conclusion

This study identifies GNLY as a potential contributor to primary open-angle glaucoma (POAG) through multi-omics analysis and experimental validation, suggesting its role in immune modulation and apoptosis. Further studies are needed to elucidate its precise mechanisms and therapeutic potential.