Pathway Analysis Integrating Genome-Wide and Functional Data Identifies PLCG2 as a Candidate Gene for Age-Related Macular Degeneration.

Purpose: Age-related macular degeneration (AMD) is the worldwide leading cause of blindness among the elderly. Although genome-wide association studies (GWAS) have identified AMD risk variants, their roles in disease etiology are not well-characterized, and they only explain a portion of AMD heritability.
Methods: We performed pathway analyses using summary statistics from the International AMD Genomics Consortium's 2016 GWAS and multiple pathway databases to identify biological pathways wherein genetic association signals for AMD may be aggregating. We determined which genes contributed most to significant pathway signals across the databases. We characterized these genes by constructing protein-protein interaction networks and performing motif analysis.
Results: We determined that eight genes (C2, C3, LIPC, MICA, NOTCH4, PLCG2, PPARA, and RAD51B) "drive" the statistical signals observed across pathways curated in the Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, and Gene Ontology (GO) databases. We further refined our definition of statistical driver gene to identify PLCG2 as a candidate gene for AMD due to its significant gene-level signals (P < 0.0001) across KEGG, Reactome, GO, and NetPath pathways. Conclusions: We performed pathway analyses on the largest available collection of advanced AMD cases and controls in the world. Eight genes strongly contributed to significant pathways from the three larger databases, and one gene (PLCG2) was central to significant pathways from all four databases. This is, to our knowledge, the first study to identify PLCG2 as a candidate gene for AMD based solely on genetic burden. Our findings reinforce the utility of integrating in silico genetic and biological pathway data to investigate the genetic architecture of AMD.

Vision loss is one of the most feared medical conditions because of its profound effect on day-to-day quality of life.1,2 Age-related macular degeneration (AMD) is the most common cause of blindness in individuals over age 60 and is responsible for almost 10% of all cases of blindness in the world.3 AMD is a late-onset disease that results from the accumulation of drusen, inflammation, and photoreceptor loss in the macular region of the eye.3 This progressive disease is categorized as either early/intermediate or advanced AMD; the latter is further subclassified as geographic atrophy (dry AMD [GA]) or choroidal neovascularization (wet AMD [CNV]).3 Early AMD is often asymptomatic and dry AMD is initially asymptomatic, but as the disease progresses, patients' central vision begins to blur and diminish.3 Wet AMD is characterized by the growth of abnormal blood vessels in the macula, which ultimately results in severe vision loss.3

Although both genetic and environmental factors shape AMD susceptibility, between 46% and 71% of the phenotypic variance of the disease is attributable to genetic factors.4 To understand the genetic architecture of AMD, the International Age-Related Macular Degeneration Genomics Consortium (IAMDGC) performed a large-scale genome-wide association study (GWAS) for advanced AMD cases and controls. They identified 52 independent genetic variants across 34 susceptibility loci for advanced AMD that are estimated to explain nearly two thirds of AMD heritability.5 Therefore, about one third of AMD heritability is still unexplained by the known loci. Although other studies have identified additional risk loci with modest effect for advanced AMD,6,7 more comprehensive approaches beyond GWAS must be used to find the remaining heritable variation for AMD.

Rather than investigating associations between single genetic variants and a phenotype, pathway analysis of GWAS data interrogates alterations in biological pathways for a trait of interest. Generally, this is done by aggregating summary statistics for these variants into genes, which are then grouped into pathways based on data in curated pathway databases.8 We hypothesize that applying this more comprehensive approach may help elucidate the genetic etiology of advanced AMD that has been indiscernible from GWAS. In this study, we performed in silico pathway analysis using the Pathway Analysis by Randomization Incorporating Structure (PARIS) software to identify biological pathways and processes enriched in genetic variation potentially associated with AMD in individuals of European descent. Because nomenclature, foci, and definitions vary across pathway databases,9 we utilized multiple databases to complement and validate our findings. Additionally, we sought to determine the central causal genes that “drive” the statistical signals observed for significant pathways identified by PARIS.

Methods

Study Subjects and GWAS Summary Statistics

The participants for this study were previously ascertained by cohorts in the IAMDGC as described.5 This included 16,144 individuals with advanced AMD and 17,832 unaffected individuals. Of the advanced AMD cases, 3235 individuals have GA only and 10,749 have CNV only. The remaining cases have both GA and CNV. All of the cases and controls used for our analyses were of European ancestry. All participants provided informed consent, and the study protocol was approved by institutional review boards as previously described.5 Data were previously collected in accordance with the tenets of the Declaration of Helsinki. The summary statistics we analyzed in this study were obtained in the 2016 GWAS performed by the IAMDGC.5 Specifically, these data include P values for 445,115 directly genotyped common and rare variants from the advanced AMD case-control results. The genotypes for these variants were generated from an array (HumanCoreExome; Illumina, San Diego, CA, USA) that was designed with additional genome-wide and custom content for AMD.5

PARIS: Knowledge-Driven Pathway Analysis of GWAS Data

To identify biological pathways enriched in genetic variants possibly contributing to advanced AMD risk, we performed in silico pathway analysis using the PARIS v2.4 software.10,11 PARIS uses variant summary statistics from GWAS, clusters them into features defined by the linkage disequilibrium (LD) structure of the genome based on a reference catalog of common genetic variants, and assigns significance to pathways based on permutation of the genome.10,11 In our analyses, we performed 100,000 permutations. PARIS also assigns empirical P values to the genes composing a pathway based on permutation testing of features within each of the genes.10,11

We performed PARIS using multiple pathway databases, including Kyoto Encyclopedia of Genes and Genomes (KEGG),12 Reactome,13 Gene Ontology (GO),14 and NetPath.15 KEGG, Reactome, and GO databases are extensive, curated biological pathway data repositories. NetPath is a specialized database that covers signaling pathways. Pathways with a P value less than 0.0001 were prioritized for further investigation. This permutation P value was calculated using the following equation: P = (1 + b)/(1 + M), where M = the number of permutations and b is the number of randomly sampled permutation scores that are greater than the observed score. To determine if the pathway associations we observed were driven by known AMD loci, we reperformed our pathway analyses excluding variants from the 34 susceptibility loci identified by the IAMDGC (defined by the 52 genomic variants) and their proxies (r2 ≥ 0.5) within 500 kb.5

Identification of Statistical Pathway Driver Genes

Due to disparate nomenclature and composition of pathways in the databases, we identified genes that overlapped across significant pathways within a database and across databases (regardless of pathway). This served to internally validate and complement our results. To interrogate the significant signals obtained from the pathways identified by PARIS, we queried which significant (P < 0.0001) genes overlapped among the significant (P < 0.0001) pathways within a pathway database. These genes were compared across the analyses done with each of the pathway databases (KEGG, Reactome, GO, and NetPath) to find statistical driver genes that had significant signals across three or more databases for the advanced AMD results.

Protein-Protein Interaction (PPI) Network for Statistical Pathway Driver Genes

We searched the Search Tool for Recurring Instances of Neighbouring Genes (STRING) database16 version 10.5 for PPIs involving the proteins encoded by the genes identified as statistical driver genes. The STRING database is composed of known and predicted PPIs based on data from curated interactions databases, high-throughput lab experiments, coexpression, and text mining in the literature. We used the high confidence (0.700) minimum required interaction score to construct the protein-protein networks of interactions based on experimental data, database entries, and coexpression.

Motif Analysis for Statistical Pathway Driver Genes

We extracted reference genome sequences for the statistical driver genes using the UCSC Genome Table Browser.17 We included 600 nucleotides upstream from the first exon and the 5′ untranslated region (UTR) in the sequences for each gene. To identify potential sequence motifs for each of these gene sets, we utilized the Multiple Expectation Maximization (EM) for Motif Elucidation (MEME) software suite.18 Sequences were considered motifs if their lengths were between 6 and 50 nucleotides. MEME was not required to find a motif in every sequence, but motifs were required to have an E-value of 0.0001. Each motif from the gene sets was then investigated in Tomtom, which looks for transcription factors (TFs) that are associated with the motif. TF binding motifs were evaluated based on the known human TF database from JASPAR19 using HOCOMOCO.20 To validate the motifs found and to test the null hypothesis of random motifs found unrelated to the statistical driver genes, 10 permutations were run on a random gene set generator for eight genes and performed the same analyses via MEME and Tomtom. We removed motifs and TFs that appeared in both the random and actual gene sets from further analysis.

Results

In Silico Pathway Analysis

We identified several biological pathways and processes from KEGG, Reactome, GO, and NetPath databases (Table 1; Supplementary Tables S1–S4) to be significantly associated with advanced AMD using PARIS. A pathway was considered significant if it had a pathway-level P value less than 0.0001. The vast majority of pathways in the four databases were not significant (Table 1). When we reperformed our pathway analyses excluding the 34 known AMD loci,5 ∼40% of the previously significant KEGG (n = 10) and GO (n = 53) pathways and over 60% of the Reactome (n = 32) pathways remained significant (Supplementary Tables S1–S3). The single NetPath pathway that was significant in our initial analysis (Wnt; Supplementary Table S4) was no longer significant in this sensitivity analysis (P = 0.00215).

Table 1

Significantly Associated Pathways Across Multiple Pathway Databases for Advanced AMD

Database	Count of Significant Pathways	Total Entries in Database	Proportion of Significant Pathways in Database
NetPath	1	26	0.038
KEGG	25	293	0.085
Reactome	50	1,748	0.029
GO	145	12,765	0.011

Pathways were considered significant if they obtained an empirical P < 0.0001.

Statistical Driver Genes Among Advanced AMD-Associated Pathways

Because pathway structure and terminology vary across databases, we determined which genes were significantly contributing to the overall pathway signals detected by PARIS. We compared the significant genes in significant pathways from KEGG, Reactome, and GO (Fig. 1; Table 2) and identified eight such genes. Upon removing variants from our analyses that fell within the 34 known AMD susceptibility loci as defined in Supplementary Table S5 in the IAMDGC GWAS,5 we found that two genes (PPARA and PLCG2) remained statistical driver genes across associated pathways from KEGG, Reactome, and GO.

Figure 1

Comparison of significant genes from AMD-associated KEGG, Reactome, and GO pathways identified by PARIS. Eight genes demonstrated significant signals across all three comparisons and are summarized in Table 2.

Table 2

Eight Statistical Pathway Driver Genes From Significant KEGG, Reactome, and GO Pathways

Gene	Chromosome	Full Gene Name (HGNC)
Statistical pathway driver genes implicated in the 2016 IAMDGC GWAS Loci
C2	6	Complement C2
MICA	6	MHC class I polypeptide-related sequence A
NOTCH4	6	Notch receptor 4
RAD51B	14	RAD51 paralog B
LIPC	15	Lipase C, hepatic type
C3	19	Complement C3
Novel genes identified with pathway analysis with PARIS
PLCG2	16	Phospholipase C gamma 2
PPARA	22	Peroxisome proliferator activated receptor alpha

The cross-database comparison of significant genes from significantly associated pathways.

To identify evidence of PPI for the proteins encoded by the eight statistical driver genes in our analyses (C2, C3, LIPC, MICA, NOTCH4, PPARA, PLCG2, and RAD51B), we queried the STRING database. Each of these proteins have multiple binding partners identified through functional studies or in silico predictions (Fig. 2). When considering no more than 50 interaction partners for each of the eight proteins, we found three distinct clusters of PPIs (Fig. 2). One cluster connects MICA, PLCG2, LIPC, C2, C3, and other immune-related proteins (Fig. 2A); another connects NOTCH4, PPARA, and other signaling proteins (Fig. 2B); and the third contains RAD51B and other DNA repair proteins (Fig. 2C).

Figure 2

PPI network generated for the proteins encoded by the eight statistical driver genes. No more than 50 interactions from the STRING database were displayed for each input protein. This threshold of interactions enabled the connection of all eight queried proteins to a network. Three distinct networks were defined by the proteins encoded by the statistical driver genes: (A) network connecting MICA, PLCG2, LIPC, C2, C3, and other immune-related proteins; (B) network connecting NOTCH4, PPARA, and other signaling proteins; (C) network connecting RAD51B and other DNA repair proteins. Types of interaction sources include coexpression (black), experimental data (magenta), and curation in databases (cyan).

Using the MEME software suite, we identified sequence motifs with known TF binding sites near the eight statistical driver gene seqeuences from the UCSC Genome Table Browser.17 Five motifs were present for most of the statistical driver genes and contain binding sites for TFs (Table 3). Only one sequence motif ([GCA][AC][CT]AG[AT]G[CA][TGA]A[AG][AT][CA]T[CA][CG][GA]T[CG][TG][CA]A[AG]AAA[ATG][AG]AAA[AT][CA][AC]A[AC]A[AC][AT][AT]A) was near all eight statistical driver genes and contained binding sites for 12 TFs.

Table 3

Sequence Motifs With TF Binding Sites Near Statistical Driver Genes

Motif Consensus Sequence	TF	P Value	Statistical Driver Genes
G[CG][TG]TG[AT]ACC[CAT][AG]G[GT][AG]GG[CT][GT][GT][AT][GA] [CG]TT[GC]C[AT]G[TA]GAGCC[GT]AGA[TA]C[GA][CG][GT][CT] C[AT][CG]	KLF5	0.0095	C2LIPCMICANOTCH4PPARARAD51B
	KLF12	0.011
	THA11	0.012
	ZN563	0.013
	IRF2	0.013
	NFIA	0.017
	ZN449	0.019
	ELF2	0.020
	ZBTB6	0.021
	RARG	0.024
[CT][TA]G[GT]C[TC]AA[CA][AG][CT][AG][GC][TA][GC]AAACCC[CA] [GC][TA][CA][TA]C[TC]A[CT][TC][AC]AA[AG]ATA[CT][AT][AG] [AC]AAA[AT]TA[GT][TCG]	PIT1	0.0052	C2LIPCMICANOTCH4PPARARAD51B
	SOX5	0.0093
	AIRE	0.010
	CEBPE	0.011
[GA][CG]CTG[CT][AT][GA][TA]CC[CA]AGCT[AGC][CT][TA][CGA][GT] [GT][GT][AT][GC]G[CTC][TG][GA]AG[GT]CAG[GA][AT]G[AC][AC] [TGC]	MAFB	0.0089	C2LIPCMICANOTCH4PPARARAD51B
	MAFF	0.010
	HTF4	0.011
	MAFK	0.012
	FOXA2	0.012
	TFE2	0.014
	BACH2	0.021
[GA]C[CT]T[CT][GC][GA]CC[TC]CCCAAA[GC][TC]GCTGGGAT[TC] AC[AG]GGCGT[GC]A[GA]CC	TFAP4	0.0047	C2LIPCMICANOTCH4PPARAPLCG2RAD51B
	ZN322	0.0062
	ZNF41	0.011
	CRX	0.013
	ZIC3	0.015
	NKX21	0.020
	GLI3	0.024
[GCA][AC][CT]AG[AT]G[CA][TGA]A[AG][AT][CA]T[CA][CG][GA]T[CG] [TG][CA]A[AG]AAA[ATG][AG]AAA[AT][CA][AC]A[AC]A[AC][AT][AT]A	HEN1	0.0025	C2C3LIPCMICANOTCH4PPARAPLCG2RAD51B
	ZSC31	0.0029
	PKNX1	0.0034
	NKX21	0.0037
	PBX3	0.0066
	TYY1	0.010
	NR2C1	0.011
	VDR	0.014
	CREB1	0.016
	RFX2	0.021
	ATF1	0.021
	CEBPE	0.022

For each motif, we identified TFs associated with the motif sequence using Tomtom. The P value represents the strength of the match between the sequence motif identified adjacent to the statistical driver genes and the curated sequences of the TF binding motifs in the HOCOMOCO database.

We further restricted our definition of statistical pathway driver gene to include genes that also strongly contributed to AMD-associated pathways from NetPath. This enabled us to further support PLCG2 as a candidate gene for advanced AMD (Fig. 3). This gene encodes a phosphodiesterase that is involved in phosphatidylinositol signaling and several other immune, metabolic, and signaling pathways curated in KEGG, Reactome, GO, and NetPath (Fig. 3). We interrogated potential interaction partners for the PLCG2 protein by constructing a PPI network for PLCG2 using the STRING database (Fig. 4). We also determined if PLCG2 harbored any suggestive associations with AMD in the IAMDGC data. None of the P values for the 65 individual PLCG2 variants we analyzed with PARIS reach genome-wide significance (P < 5 × 10−8), but several of them (n = 14) were nominally associated (P < 0.05) with advanced AMD (Fig. 5). The single-variant association results from PLCG2 are not highly correlated based on LD structure using the 1000 Genomes Project (Fig. 5), which indicates that the concentration of nominally significant results in this gene is not merely due to LD.

Figure 3

Identification of PLCG2 as a candidate gene for advanced AMD. A comparison of the significant genes from significant KEGG, GO, NetPath, and Reactome pathways in our PARIS pathway analysis converged on one gene (PLCG2), which encodes a protein that is common to several pathways.

Figure 4

PPI network generated for PLCG2. No more than 10 interactions were displayed. Types of interaction sources include coexpression (black), experimental data (magenta), and curation in databases (cyan).

Figure 5

LocusZoom Plot of P values for the 65 PLCG2 variants in the IAMDGC advanced AMD case-control analysis. These variants were either within the gene boundaries (human genome build 37) of PLCG2 or within 50 kb of these boundaries. P values were generated by the IAMDGC in their advanced AMD case-control GWAS published in 2016.5 LD estimates (r2) are based on the European (EUR) population from the 1000 Genomes Project (November 2014 release).

Discussion

Using knowledge-driven pathway analysis on GWAS data, we uncovered pathways that were enriched in variation potentially associated with AMD in individuals of European descent. Our study is, to our knowledge, the first to perform such analyses on the largest available advanced AMD case-control association dataset. We found several signaling, immune, metabolic, and disease-related pathways from the KEGG, Reactome, GO, and NetPath databases that are associated with advanced AMD. Our sensitivity analysis demonstrated that several of the pathways from KEGG, Reactome, and GO (Supplementary Tables S1–S3) remained associated with advanced AMD following the exclusion of the 34 AMD susceptibility loci described earlier.5 This suggests that modest effects aggregating in these pathways may contribute to the missing heritability of AMD. Although the Wnt pathway from NetPath was no longer significant in our sensitivity analysis, the Wnt signaling pathway from GO remained associated with AMD. This results from the difference in the pathway definitions. These pathways are nearly identical in size (n = 45 and 41 genes for NetPath and GO, respectively); however, only two genes overlap between them (PLCG2 and FZD4). Furthermore, the Wnt signaling pathway in KEGG (n = 140 genes) and the signaling by Wnt pathway in Reactome (n = 294 genes) only achieved pathway-level P values of 0.032 and 0.037 in our analyses, respectively. These pathway definition differences further justify our use of multiple curated databases in our analyses to uncover AMD-associated pathways and genes driving their statistical significance.

Due to varying nomenclature for pathways across databases and as a way of internal validation, we focused on eight statistical driver genes (C2, C3, LIPC, MICA, NOTCH4, PPARA, PLCG2, and RAD51B) that were consistently significant across GO, Reactome, and KEGG pathways. PPARA and PLCG2 were not previously identified as a part of the 34 IAMDGC loci associated with AMD risk. The strongest single-marker P values observed in PLCG2 and PPARA were 2.05 × 10−4 and 3.10 × 10−5, respectively, and do not meet the classical GWAS significance levels. In our sensitivity analysis, PPARA and PLCG2 remained statistical driver genes in pathways from KEGG, Reactome, and GO, suggesting that pathway analysis can identify novel AMD genes. Additionally, the aggregation of nominally significant independent variants in PLCG2 suggests that the gene-wide significance of PLCG2 is greater than that of the individual variants and emphasizes the power of pathway analysis for identifying gene-wide signals rather than single-variant associations.

DNA motif analysis identified five sequence motifs adjacent to the eight statistical driver genes in their promoter regions. These motifs represent sites of known TF binding and suggest that the expression of these genes may be controlled by similar mechanisms. One motif ([GCA][AC][CT]AG[AT]G[CA][TGA]A[AG][AT][CA]T[CA][CG][GA]T[CG][TG][CA]A[AG]AAA[ATG][AG]AAA[AT][CA][AC]A[AC]A[AC][AT][AT]A) was adjacent to the start positions of all eight statistical driver genes and contains known binding sites of several TFs (Table 3). Functional studies are required to confirm these in silico findings and elucidate the transcriptional mechanisms of these statistical driver genes in the context of AMD.

One gene, PLCG2, was central to multiple pathways in all four databases and remained significant after our sensitivity analysis. PLCG2 encodes a signaling enzyme (phospholipase C gamma 2, PLCG2) that utilizes calcium to catalyze the hydrolysis of PIP2 into second messengers IP3 and DAG.21 These molecules initiate intracellular calcium flux and activate protein kinase C, respectively.21 The enzymatic activity of PLCG2 results from tyrosine phosphorylation performed by growth factor receptors, immune receptors, and G protein-coupled receptors as well as the activity of lipid-derived second messengers in the cell.21 This enzyme is highly expressed in cells of hematopoietic origin and is responsible for regulating immune responses and platelet adhesion and spreading.22–26

The PLCG2 protein interacts with several members (HCK, LYN, PIK3R1, and SYK) of the microglia pathogen phagocytosis pathway in humans.27 Its interaction partners also play roles in oxidative stress, angiogenesis, and platelet activation. BLNK and BTK are central to facilitating B-cell apoptosis following oxidative stress.28,29 Exposure to oxidative stress activates EGFR, which promotes retinal epithelial cell health and survival through EGFR/Akt, PI3K, and ERK/MAPK signaling pathways.30,31 EGFR downstream signaling also contributes to retinal pigment epithelial cell proliferation and migration in wound healing.32,33 PIK3R1 is a regulatory subunit of PI3K in the PI3K/Akt/mTOR pathway, which is a possible target for treating ocular neovascularization.34 PI3K and Tec protein kinases regulate platelet activation,35 and signaling cascades from LCP2 (also called SLP-76) and SYK are responsible for separating blood and lymphatic vasculatures in the human body.36 These interactions and processes, coupled with PLCG2's role in the VEGF pathway,37,38 could be pertinent for understanding the role of PLCG2 and its interaction partners in the choroidal neovascularization subtype of advanced AMD. In the CNV-only case-control GWAS performed by the IAMDGC, no PLCG2 variants were genome-wide significant; however, 13 variants were nominally associated with CNV (P < 0.05).5 Of the 65 PLCG2 variants analyzed by PARIS, 31 exhibited lower P values in the CNV-specific IAMDGC GWAS than in the combined advanced AMD IAMDGC GWAS.

Heterozygous gain-of-function mutations in PLCG2 result in constitutive phospholipase activity and PLCG2-associated antibody deficiency and immune dysregulation, which is characterized by immunodeficiency and autoimmunity.39 This gene was recently identified as a candidate gene for rheumatoid arthritis (RA) due to its overexpression in RA patients compared to controls.40 Genetic risk scores for RA are associated with increased AMD risk,41 and individuals with RA are at a higher risk of developing AMD.42 PLCG2 is also highly expressed in microglia43 and has been previously implicated in the genetic etiology of late-onset Alzheimer's disease (LOAD).44,45 Specifically, GWAS identified a protective effect for a rare variant in the coding region of PLCG2 on LOAD.44,45 This variant is considered hypermorphic because the mutant enzyme experiences a small increase in enzymatic activity compared to wild-type enzyme, which would imply that mildly activating PLCG2 could be a therapeutic intervention for LOAD.43 Functional studies would need to be performed to determine if PLCG2's enzymatic activity could be modulated by a similar mechanism in patients with AMD.

Although PLCG2 has not been previously associated with AMD in a case-control GWAS, variants in this gene were associated with AMD when accounting for birth control pill usage in women with CNV.46 These associations were undetectable when gene-environment interactions between PLCG2 variants and exogenous estrogen exposure were not considered.46 Other interaction studies have identified PLCG2 variants as genetic modifiers of previously identified associations among menopausal hormone therapy, mammographic density, and breast cancer risk, which could suggest sex-specific effects of genetic variants in this gene for disease risk.47,48

While our study provides in silico evidence for the roles of these statistical driver genes and pathways in AMD, it does not biologically confirm them. Functional studies are required to determine causality for these genes and pathways in patients with AMD. Knowledge-driven pathway analyses are subject to the quality and coverage of the knowledge in a given database. We attempted to circumvent this limitation by utilizing multiple databases in our analyses and integrating our results. The GWAS data used in this study were generated from individuals of European descent. Consequently, these findings may not be applicable to non-European populations. The IAMDGC GWAS dataset is considered the largest available dataset for advanced AMD cases and controls in the world. We are unaware of any comparable datasets available for replication.

Click here for additional data file.