Mendelian Randomization and Transcriptome-Wide Association Analysis Identified Genes That Were Pleiotropically Associated with Intraocular Pressure

Background: Intraocular pressure (IOP) is a major modifiable risk factor for glaucoma. However, the mechanisms underlying the controlling of IOP remain to be elucidated. Objective: To prioritize genes that are pleiotropically associated with IOP. Methods: We adopted a two-sample Mendelian randomization method, named summary-based Mendelian randomization (SMR), to examine the pleiotropic effect of gene expression on IOP. The SMR analyses were based on summarized data from a genome-wide association study (GWAS) on IOP. We conducted separate SMR analyses using Genotype-Tissue Expression (GTEx) and Consortium for the Architecture of Gene Expression (CAGE) expression quantitative trait loci (eQTL) data. Additionally, we performed a transcriptome-wide association study (TWAS) to identify genes whose cis-regulated expression levels were associated with IOP. Results: We identified 19 and 25 genes showing pleiotropic association with IOP using the GTEx and CAGE eQTL data, respectively. RP11-259G18.3 (PSMR = 2.66 × 10−6), KANSL1-AS1 (PSMR = 2.78 × 10−6), and RP11-259G18.2 (PSMR = 2.91 × 10−6) were the top three genes using the GTEx eQTL data. LRRC37A4 (PSMR = 1.19 × 10−5), MGC57346 (PSMR = 1.19 × 10−5), and RNF167 (PSMR = 1.53 × 10−5) were the top three genes using the CAGE eQTL data. Most of the identified genes were found in or near the 17q21.31 genomic region. Additionally, our TWAS analysis identified 18 significant genes whose expression was associated with IOP. Of these, 12 and 4 were also identified by the SMR analysis using the GTEx and CAGE eQTL data, respectively. Conclusions: Our findings suggest that the 17q21.31 genomic region may play a critical role in the regulation of IOP.


Introduction
Intraocular pressure (IOP) refers to the fluid pressure inside the eye, which is affected by the production and drainage of the aqueous humor. Excessive production or insufficient drainage of aqueous humor can lead to ocular hypertension (OHT), which has a prevalence of approximately 1.6% in the general population over 30 years old [1]. The prevalence of OHT can be higher in persons over 40 years old, ranging from 2.6% to 3.6% [2][3][4]. IOP can be modified by various treatments, such as eye drops, laser, and surgery [5]. Previous studies have shown that each 1 mm Hg drop of IOP can lower the risk of glaucoma progression by 10% [6,7]. Therefore, IOP represents a major modifiable risk factor of glaucoma [8], the most common cause of irreversible blindness worldwide [9]. glaucoma progression by 10% [6,7]. Therefore, IOP represents a major modifiable risk factor of glaucoma [8], the most common cause of irreversible blindness worldwide [9].
In this paper, we investigated genes that are pleiotropically associated with IOP by using a recently developed two-sample Mendelian randomization (MR), named summary-based Mendelian randomization (SMR), that integrates genome-wide association study (GWAS) summary data for IOP and cis-eQTL (expression quantitative trait loci) data. We also conducted a transcriptome-wide association study (TWAS) to identify genes whose cis-regulated expression levels are associated with IOP.

Editorial Policies and Ethical Considerations
This study utilized GWAS summary results for IOP and eQTL data. All the data were publicly available. As a result, ethical considerations are not needed. The analytical process of the present study is illustrated in Figure 1.

GWAS Data for IOP
The GWAS summary data for IOP were provided by a recent multi-trait genomewide association meta-analysis of optic disc parameters [28]. The results were based on genetic data imputed using the Haplotype Reference Consortium (HRC) and included 11 cohorts from the International Glaucoma Genetics Consortium (IGGC). A total of 31,269 participants of European ancestry were included in the meta-analysis. An additive genetic model was assumed by all the participating studies, adjusting for age, sex, and at least the first five principal components, as well as cohort-specific covariates when necessary. The GWAS Figure 1. Flow chart for the SMR analyses. (A) SMR analysis using GTEx eQTL data for blood; and (B) SMR analysis using CAGE eQTL data for blood. CAGE, Consortium for the Architecture of Gene Expression; eQTL, expression quantitative trait loci; GWAS, genome-wide association studies; GTEx, Genotype-Tissue Expression; LD, linkage disequilibrium; IOP, intraocular pressure; SMR, summary-based Mendelian randomization; SNP, single nucleotide polymorphisms.

GWAS Data for IOP
The GWAS summary data for IOP were provided by a recent multi-trait genomewide association meta-analysis of optic disc parameters [28]. The results were based on genetic data imputed using the Haplotype Reference Consortium (HRC) and included 11 cohorts from the International Glaucoma Genetics Consortium (IGGC). A total of 31,269 participants of European ancestry were included in the meta-analysis. An additive genetic model was assumed by all the participating studies, adjusting for age, sex, and at least the first five principal components, as well as cohort-specific covariates when necessary. The GWAS summarized data are publicly available and can be downloaded from http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST009001-GCST0 10000/GCST009413 (accessed on 26 July 2021).

eQTL Data
The SMR analyses utilized cis-eQTL genetic variants as instrumental variables (IVs) for gene expression. Because eQTL data for the eye were unavailable, eQTL data for blood from Genotype-Tissue Expression (GTEx) and Consortium for the Architecture of Gene Expression (CAGE) were used for the SMR analyses. The V7 release of the GTEx eQTL data for blood [29] included 338 participants, while the CAGE eQTL data for blood [30] included 2765 participants. The eQTL data can be downloaded at https://cnsgenomics.com/data/ SMR/#eQTLsummarydata (accessed on 26 July 2021).

SMR Analysis
The Mendelian analyses were conducted using the method implemented in the SMR software version 1.3.1 [31]. SMR applies the principles of MR, integrating GWAS and eQTL summary statistics to explore the pleiotropic association between gene expression and a trait. In SMR, genetic variants are used as the IVs, and estimation of the pleiotropic association can be made because the inherited genetic variants are independent of potentially confounding factors [32]. Here, 'pleiotropic association' between gene expression and IOP refers to either pleiotropy (i.e., both gene expression and IOP are affected by the same causal variant) or causality (i.e., the effect of a causal variant on IOP is mediated by gene expression). In either case, gene expression is the exposure but not an outcome (i.e., SNP > gene expression). SMR is not designed to examine whether a trait influences gene expression.
The SMR analysis followed a similar approach as described in a previous publication [33], using all the default settings in SMR software. Details about the settings adopted in the SMR analyses can be found in Table S1. The multiple testing was adjusted by the false discovery rate (FDR).

TWAS Analysis
We also conducted a TWAS analysis [34] integrating the GWAS summary statistics of IOP and pre-computed gene expression weights to further explore genes whose cisregulated expression is associated with IOP. This approach imputed gene expression from summarized GWAS data to test its association with the phenotype of interest. The effect size of the expression of a gene for a trait can be viewed as a linear combination of genetic effects on a trait, with weights calculated based on the correlation between SNPs and gene expression while accounting for linkage disequilibrium (LD) among the SNPs. The weights are generally pre-computed using data from a relatively small set of reference individuals for whom both gene expression and genetic variations (SNPs) are available [34]. Although this TWAS approach is conceptually similar to SMR, it is less strict than SMR in that it aims to identify genes whose genetically controlled expression is associated with a disease, while SMR aims to identify genes that are pleiotropically associated with a disease. As we observed multiple significant genes located in 17q21.31, for this risk locus, we performed joint/conditional analysis, a process to identify which genes represent independent associations (i.e., jointly significant), and which genes are not significant after accounting for the predicted expression of other genes in the region (i.e., marginally significant). The analyses were performed by using FUSION software. We applied the weights that were pre-computed from the GTEx v7 whole-blood reference expression panel and adopted all the default settings in FUSION.

Basic Information of the Summarized Data
The CAGE eQTL data have a much larger number of participants compared to the GTEx eQTL data (2765 vs. 338), as well as a larger number of eligible probes (8524 vs. 4543). After checking the allele frequencies between the datasets and performing LD pruning, there were more than 6.2 million eligible SNPs in each SMR analysis (Table 1).

Pleiotropic Association with IOP
Using the GTEx eQTL data, our SMR analysis identified 19 genes that are pleiotropically associated with IOP (Tables 2 and S2 The GWAS summarized data were provided by the study of Bonnemaijer et al. [28] and can be downloaded at http://ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST009001-GCST010000/GCST009413 (accessed on 26 July 2021). The CAGE and GTEx eQTL data can be downloaded at https://cnsgenomics.com/data/SMR/#eQTLsummarydata (accessed on 26 July 2021). Only the top ten significant results are shown for each SMR analysis. For the complete significant results, please see Tables S1 and S2. Top SNP means the top cis-eQTL and was used as the instrument variable. PeQTL is the p-value of the top associated cis-eQTL in the eQTL analysis, and PGWAS is the p-value for the top associated cis-eQTL in the GWAS analysis. Beta is the estimated effect size in the SMR analysis, SE is the corresponding standard error, PSMR is the p-value for the SMR analysis, PHEIDI is the p-value for the HEIDI test, and the Q value is the adjusted p-value found using FDR. The FDR was calculated at p = 10 −3 threshold. CHR, chromosome; HEIDI, heterogeneity in dependent instruments; IOP, intraocular pressure; SNP, single-nucleotide polymorphism; SMR, summary-based MendelianrRandomization; QTL, quantitative trait loci; FDR, false discovery rate; GWAS, genomewide association studies.    The GWAS summarized data were provided by the study of Bonnemaijer et al. [28] and can be downloaded at http: //ftp.ebi.ac.uk/pub/databases/gwas/summary_statistics/GCST009001-GCST010000/GCST009413 (accessed on 26 July 2021). The CAGE and GTEx eQTL data can be downloaded at https://cnsgenomics.com/data/SMR/ #eQTLsummarydata (accessed on 26 July 2021). Only the top ten significant results are shown for each SMR analysis. For the complete significant results, please see Tables S1 and S2. Top SNP means the top cis-eQTL and was used as the instrument variable. P eQTL is the p-value of the top associated cis-eQTL in the eQTL analysis, and P GWAS is the p-value for the top associated cis-eQTL in the GWAS analysis. Beta is the estimated effect size in the SMR analysis, SE is the corresponding standard error, P SMR is the p-value for the SMR analysis, P HEIDI is the p-value for the HEIDI test, and the Q value is the adjusted p-value found using FDR. The FDR was calculated at p = 10 −3 threshold. CHR, chromosome; HEIDI, heterogeneity in dependent instruments; IOP, intraocular pressure; SNP, single-nucleotide polymorphism; SMR, summary-based MendelianrRandomization; QTL, quantitative trait loci; FDR, false discovery rate; GWAS, genome-wide association studies.

Discussion
In this study, we identified multiple genes showing a pleiotropic association with IOP through SMR and TWAS approaches. Our results confirmed findings from previous studies and revealed novel genes related to IOP regulation.
In our SMR analyses, the IVs were based on eQTL data, and the exposure was (transcriptome-wide) gene expression. The GWAS used genetic data imputed based on the HRC and included 11 cohorts from the International Glaucoma Genetics Consortium (IGGC) [28]. The GTEx eQTL data were based on deceased donors [29], and the CAGE eQTL data were based on samples from five cohorts [30]. These cohorts were not part of the GWAS summary results for IOP. Therefore, there is no overlap between the samples. MR can be conducted based on data from one sample or two samples. The summary association results came from the same individuals in the one-sample MR, and from different, potentially overlapping sets of individuals in the two-sample MR [35]. We chose the Figure 6. Joint/conditional analysis of TWAS significant loci on 17q21.31. The top panel of the joint/conditional plot displays all genes that are in the loci (usually gray), with marginally significant TWAS association genes highlighted in blue, and jointly significant genes in green. The bottom panel is the Manhattan plot of the original GWAS summary statistics data before (gray) and after (blue) conditioning on the green genes. TWAS, transcriptome-wide association study; GWAS, genome-wide association study.

Discussion
In this study, we identified multiple genes showing a pleiotropic association with IOP through SMR and TWAS approaches. Our results confirmed findings from previous studies and revealed novel genes related to IOP regulation.
In our SMR analyses, the IVs were based on eQTL data, and the exposure was (transcriptome-wide) gene expression. The GWAS used genetic data imputed based on the HRC and included 11 cohorts from the International Glaucoma Genetics Consortium (IGGC) [28]. The GTEx eQTL data were based on deceased donors [29], and the CAGE eQTL data were based on samples from five cohorts [30]. These cohorts were not part of the GWAS summary results for IOP. Therefore, there is no overlap between the samples. MR can be conducted based on data from one sample or two samples. The summary association results came from the same individuals in the one-sample MR, and from different, potentially overlapping sets of individuals in the two-sample MR [35]. We chose the 'two-sample MR' over the 'one-sample MR' for several reasons: (1) The eQTL data are unavailable for the subjects in the GWAS data; (2) using the association results from the same or partially overlapping samples may introduce weak instrument bias [35]; and, (3) the power of an MR can be greatly increased by using a two-sample MR approach [36,37].
In our study, most of the genes showing a pleiotropic association with IOP are near 17q21.31 (chr17: 40900001-44900000, GRCh37/hg19), a structurally complex and evolutionarily dynamic region of the genome [38][39][40]. This region contains a~970 kb inversion of the MAPT locus in populations of European ancestry [41]. MAPT (microtubule-associated protein tau) is associated with both the ganglion cell inner plexiform layer (GCIPL) and the retinal nerve fiber layer (RNFL), indicating that it might impact glaucoma pathogenesis through modulation of retinal thickness [42]. The MAPT locus has two divergent haplotypes, H1 (direct orientation) and H2 (inverted orientation), with distinct functional impacts [38]. Although both GTEx and CAGE do have eQTL data for MAPT, it was dropped in the SMR analysis because SMR only includes cis-eQTL with a p-value < 5 × 10 −8 (The minimum p-value is 1.41 × 10 −4 for cis-eQTL in the GTEx and 2.62 × 10 −7 in the CAGE). Despite that, some of the genes identified in our study were reported to be associated with MAPT haplotypes. For example, LRRC37A4 is the top-hit gene in the SMR analysis using the CAGE eQTL data, and the H1 haplotype of MAPT is associated with an increased expression of it [41]. Moreover, several other identified genes in the 17q21.31 region are either associated with IOP or act collectively in influencing IOP or associated traits. For instance, a GWAS identified 139 genetic loci associated with the macular thickness (MT), including genetic variants in KANSL1, LRRC37A4P-MAPK8IP1P2, and NSF [43]. In addition, KANSL1-AS1 (identified in GTEx), LRRC37A2 (identified in CAGE), and OR7E14P were found to form a regulatory cluster in influencing both IOP and MT [44]. Together, these findings suggest that 17q21.31 is crucial for IOP regulation. However, the exact pathogenesis remains unclear. Further investigation is needed to elucidate the exact functions of this region and examine its biological role in influencing IOP and the pathogenesis of glaucoma.
We found a significant pleiotropic association between ABO and IOP using GTEx eQTL data. ABO (Alpha 1-3-N-Acetylgalactosaminyltransferase and Alpha 1-3-Galactosyltransferase) is located on chromosome 9q34.2 and encodes proteins related to the first discovered blood group system [45]. It has seven exons and six introns [46]. A variation in ABO forms the basis of the ABO blood group [47]. Genetic variants in ABO are associated with various health conditions, such as diabetes, thromboembolism, myocardial infractions, atherosclerosis, and stroke [48,49]. In a previous multi-ancestry meta-analysis, the International Glaucoma Genetics Consortium (IGGC) revealed a novel ABO polymorphism (rs8176693) associated with IOP [50]. A later meta-analysis showed that the ABO polymorphism rs8176741 was significantly associated with IOP, vertical cup-disc ratio (VCDR), and cup area [51]. Despite these encouraging findings, the exact mechanisms underlying the observed association between genetic variants in ABO and IOP remain to be elucidated. More research is needed to understand the functions and roles of ABO in influencing IOP.
In our study, we found that some genes, such as AFAP1, showed a pleiotropic association with IOP. AFAP1 (Actin Filament Associated Protein 1) is located on 4p16.1 and encodes a protein that binds and crosslinks filaments [52,53]. Actin cytoskeleton-modulating signals are involved in the regulation of aqueous outflow and IOP [54,55]. Two SNPs in AFAP1 (rs4619890 and rs11732100) have been reported to be associated with POAG in GWAS studies [56,57]. In European-ancestry populations, the two SNPs are moderately associated with another SNP in AFAP1 (rs28795989) linked to IOP [25]. Another GWAS found that rs6816389 in AFAP1 is associated with IOP in European-ancestry participants [58]. Moreover, the expression of AFAP1 was detected in the trabecular meshwork, retina (including retinal ganglion cells [RGCs]), and optic nerve of a normal human eye and a glaucoma-tous eye [57]. Together, existing evidence implies AFAP1's potential involvement in the pathogenesis of glaucoma.
We also identified several other genes that are not on 17q21.31 but show a significant pleiotropic association with IOP regulation, such as SGTB and TEF. The SGTB gene, also known as the small glutamine-rich tetratricopeptide repeat (TPR)-containing beta, is located on 5q12.3 and belongs to the SGT (small glutamine-rich TPR-containing protein) family. SGT proteins have been associated with a variety of biological processes, including neuronal synaptic transmission, cell cycle regulation, protein folding, and apoptosis [59][60][61]. Previous research discovered that SGTB interacts with Brother of CDO (BOC) and modulates its surface presence. This subsequently leads to JNK activation, which, in turn, facilitates neuronal differentiation and the growth of neurites [62]. It was found that genetic variants near SGTB are related to IOP or corneal thickness (CCT) [63,64]. The TEF gene, also known as thyrotroph embryonic factor, is located on 22q13.2. The TEF protein belongs to the PAR (proline and acidic amino acid-rich) subfamily of bZIP transcription factors [65]. The proteins encoded by these genes have recently been demonstrated to regulate the expression of many enzymes and molecules involved in detoxification and drug metabolism [8]. A previous study found that the genetic polymorphism rs6519240 near TEF is associated with refractive error [66]. However, studies on the relationship between SGTB, TEF, and IOP regulation are scarce. Further investigation is necessary to elucidate the role of these genes in IOP regulation.
An SMR analysis relies on three key assumptions. First, the genotype is associated with gene expression. Second, confounding factors that bias the associations between gene expression and IOP are not associated with the genotype. Third, the genotype is related to IOP only through its association with gene expression. For assumption 1, we used the default p-value threshold of 5 × 10 −8 to select the top associated eQTL in our SMR analyses. Therefore, the genetic variants selected as IVs are indeed strongly associated with gene expression, and weak instruments are unlikely to be a big concern. Assumption 2 is difficult to verify directly, as SMR analyses use summarized data. The assumption is often based on the biological belief that genotypes are unrelated to confounding factors, such as socioeconomics and behavioral characteristics [32]. For assumption 3, horizontal pleiotropy was found in about half of the significant causal relationships in MR, which could introduce distortions as high as 201% in the causal estimates. Horizontal pleiotropy could induce false-positive causal findings in up to 10% of the relationships [67]. We observed some pleiotropic associations with significant HEIDI tests, suggesting horizontal pleiotropy (Table 2). Therefore, caution should be exercised when interpreting the corresponding findings.
Our study has limitations. The eQTL data are based on limited sample sizes, which may affect the statistical power. Additionally, the eQTL data have limited eligible probes. As a result, we may have missed some important genes. Despite this, the power of SMR has been examined through extensive simulation studies. The simulations showed that SMR was equivalent to MR analysis if the genotype, gene expression, and phenotype data were from the same sample. The power of SMR could be greatly increased if the eQTL data and GWAS summary results were from two independent samples with very large sample sizes [31]. We believe that concerns about the power of our SMR analyses are minimal, especially in the SMR analysis using the CAGE eQTL data. The SMR approach cannot differentiate between pleiotropy and causality. We used eQTL data from blood because usable eQTL data from the eye are unavailable. Future studies using eye eQTL data are needed to validate our findings. Our SMR analyses used data from participants of European ancestry. Since the prevalence of POAG is ethnicity-specific, it is reasonable to postulate that the GWAS results might also be ethnicity-specific. Therefore, our results might not be generalized to other ethnic populations.

Conclusions
We identified several genes that are pleiotropically associated with IOP. Our results indicate that the 17q21.31 genomic region could be crucial for IOP regulation. Future studies are necessary to clarify the collective actions of the genes identified in the 17q21.31 genomic region and the roles of the identified genes on the other chromosomes in the IOP regulation.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/genes14051027/s1, Figure S1: Pleiotropic association of ABO with IOP; Figure S2: Pleiotropic association of SGTB with IOP; Table S1: Default parameters in the SMR analyses; Table S2: Genes showing significantly pleiotropic association with IOP using GTEx eQTL data; Table S3: Genes showing significantly pleiotropic association with IOP using CAGE eQTL data; Table S4: Genes that showed transcriptome-wide significant associations with IOP. Institutional Review Board Statement: Ethical review and approval were waived for this study because the study only used publicly available summary data.

Informed Consent Statement: Not applicable.
Data Availability Statement: All the data generated or analyzed during this study are publicly available as specified in Section 2 of this paper. Specifically, the eQTL data can be downloaded at https: //cnsgenomics.com/data/SMR/#eQTLsummarydata (accessed on 26 July 2021), and the GWAS summarized data can be downloaded at http://ftp.ebi.ac.uk/pub/databases/gwas/summary_ statistics/GCST009001-GCST010000/GCST009413 (accessed on 26 July 2021).

Conflicts of Interest:
The authors declare no conflict of interest.