Characterizing the Causal Pathway for Genetic Variants Associated with Neurological Phenotypes Using Human Brain-Derived Proteome Data.

Leveraging high-dimensional molecular datasets can help us develop mechanistic insight into associations between genetic variants and complex traits. In this study, we integrated human proteome data derived from brain tissue to evaluate whether targeted proteins putatively mediate the effects of genetic variants on seven neurological phenotypes (Alzheimer disease, amyotrophic lateral sclerosis, depression, insomnia, intelligence, neuroticism, and schizophrenia). Applying the principles of Mendelian randomization (MR) systematically across the genome highlighted 43 effects between genetically predicted proteins derived from the dorsolateral prefrontal cortex and these outcomes. Furthermore, genetic colocalization provided evidence that the same causal variant at 12 of these loci was responsible for variation in both protein and neurological phenotype. This included genes such as DCC, which encodes the netrin-1 receptor and has an important role in the development of the nervous system (p = 4.29 × 10-11 with neuroticism), as well as SARM1, which has been previously implicated in axonal degeneration (p = 1.76 × 10-08 with amyotrophic lateral sclerosis). We additionally conducted a phenome-wide MR study for each of these 12 genes to assess potential pleiotropic effects on 700 complex traits and diseases. Our findings suggest that genes such as SNX32, which was initially associated with increased risk of Alzheimer disease, may potentially influence other complex traits in the opposite direction. In contrast, genes such as CTSH (which was also associated with Alzheimer disease) and SARM1 may make worthwhile therapeutic targets because they did not have genetically predicted effects on any of the other phenotypes after correcting for multiple testing.

Introduction

The widespread application of genome-wide association studies (GWAS) over the last decade has drastically advanced the discovery of genetic variants associated with complex traits and disease., However, the underlying biological mechanisms responsible for the vast majority of these effects have been challenging to decipher. Alterations to circulating protein levels are likely to reside along the causal pathway between genetic variant and phenotypic variation., This has led to recent studies characterizing genetic variants associated with protein levels (known as protein quantitative trait loci [pQTL]) by integrating their effects on traits with results from GWAS.6, 7, 8, 9, 10 Furthermore, findings from these endeavors can be valuable for drug target prioritization, particularly given that therapeutic targets with support from genetic association studies are more likely to succeed in clinical trials., These estimates have recently been revised, suggesting that support from human genetics can increase approval rates for drugs targeting GWAS traits by over 2-fold.

Circulating protein levels, as with other molecular traits such as gene expression and epigenetic processes, are known to vary depending on the tissue type they are investigated in.14, 15, 16 This has therefore been a limitation for previous studies using human pQTL data, which have typically been confined to plasma proteins derived from whole blood., Although these data are more readily accessible in larger samples, they may not necessarily capture biological effects from tissue types that are more relevant to the disease being studied. For example, it would be expected that data derived from brain tissue would be the most pertinent for characterizing genetic variants associated with neurological phenotypes. This is because, due to their functionality in the brain, the underlying genes responsible for these effects are likely to exert their influence on cognitive traits and psychiatric disorders.,

In this study, we have leveraged pQTL data derived from human brain tissue through the use of the brainQTL resource and findings from large-scale GWAS of seven neurological phenotypes (Alzheimer disease, amyotrophic lateral sclerosis, depression, insomnia, intelligence, neuroticism, and schizophrenia) (Table S1).21, 22, 23, 24, 25, 26, 27 These effects were integrated using a Mendelian randomization (MR) framework which harnesses genetic colocalization to highlight loci where dorsolateral prefrontal cortex proteins and neurological phenotypes are influenced by a shared causal variant. We conducted in-depth evaluations of proteins identified in this analysis by undertaking a phenome-wide MR study to assess their association with 700 complex traits and diseases.

Material and Methods

Data Resources

All genetic effects onto brain-derived protein levels were downloaded from the brainQTL resource (see Web Resources). Further details can be found in the brainQTL paper. In brief, genotype and proteome data on 7,901 total proteins were available from 144 post-mortem samples from the Religious Orders Study (ROS) and the Memory and Aging Project (MAP). Access to full datasets from the ROS and MAP studies can be requested at the URL in the Web Resources section. pQTL identification was restricted to variants with a minor allele frequency of over 5% and which resided within 100kbs of protein coding genes based on the UCSC genome browser (build GRCh37/hg19).

We obtained genome-wide summary statistics for seven phenotypes by using findings from GWAS for which summary statistics have been made available: Alzheimer disease, amyotrophic lateral sclerosis, depression, insomnia, intelligence, neuroticism, and schizophrenia. More detailed information on all GWAS datasets can be found in Table S1.

Statistical Analysis

Linkage disequilibrium (LD) clumping was undertaken to identify independent pQTL for downstream analyses. This was achieved using PLINK based on a r2 < 0.01 with a reference panel consisting of 10,000 unrelated individuals from the UK Biobank study who were of European descent., This reference panel was selected because individuals from the ROS and MAP study were reported to be of European descent. We generated F-statistics for instruments as proposed previously by Bowden et al.:

where γj is the SNP-exposure association and σXj is the standard deviation for the SNP-exposure association for variant j. Instruments with an F-statistic > 15 were selected to reduce the likelihood of weak instrument bias in downstream analyses.

MR estimates were derived based on the Wald Ratio method using the “TwoSampleMR” package. Estimates were then filtered based on a multiple testing threshold of p < 0.05/number of proteins analyzed. Testing of seven neurological phenotypes was not taken into account in this correction, and although some of these share genetic architecture (e.g., schizophrenia and depressive symptoms have previously been reported to have an LD score regression coefficient of rG = 0.82), it was important to reinforce results surviving this cut-off with evidence of genetic colocalization. We therefore carried forward the loci that survived our correction threshold from the previous analysis, and we analyzed them with the “coloc” R package and eCAVIAR method using default parameters. Evidence of genetic colocalization was defined as having either a posterior probability of association (PPA) > 0.8 from the coloc method or a colocalization posterior probability (CLPP) > 0.01 for eCAVIAR (both of these were proposed by the developers of these methods). Results that provide evidence based on these parameters suggest that at these loci, brain-derived proteins and neurological phenotypes share a common causal variant.

We subsequently applied this analysis pipeline to 700 complex traits and diseases (a full list of which are in Table S5) by using the TwoSample MR package for any predicted effects which also had evidence of genetic colocalization from the previous analysis. These phenotypes were selected based on the following criteria: (A) outcomes which provided evidence of heritability (based on a heuristic Pheritability < 0.05) according to analyses undertaken previously by the Neale Lab (see Web Resources) or (B) outcomes from GWAS consortia based on previously defined criteria. Specifically, this included outcomes analyzed by GWAS studies that reported betas, standard error, and effect alleles for over 100,000 genetic variants and that were undertaken in a European population of over 1,000 individuals.

We also performed additional analyses for SARM1 and CTSH, for which our initial analysis did not provide evidence of any putative side effects based on Bonferroni corrections. This was undertaken by comparing the estimates from the analysis of these targets on the 700 outcomes with 500 randomly selected pQTL from our initial sample. We then ranked MR effect estimates and calculated permutated p values based on the rank of SARM1 and CTSH effects compared to those of the other 500 proteins. Any effects which provided evidence based on Ppermuted < 0.05 were further subject to the same genetic colocalization analysis as previously undertaken. Lastly, we investigated evidence of genetic colocalization between pQTL and expression quantitative trait loci (eQTL), also derived from brain tissue, based on analyses by Qi et al. (n = 1,194).

Results

A Proteome-wide Mendelian Randomization Study of Neurological Traits and Psychiatric Disorders Using Brain-Derived Data

Applying our selection critieria to identify independent pQTL resulted in 692 proteins which were eligible for analysis (Table S2). However, all of these proteins could only be instrumented using a single pQTL, and therefore all MR estimates derived in this study are based on the Wald ratio method. As sample sizes of protein data derived from disease-relevant tissues increase in the future, the MR framework proposed here has the capacity to harness multiple genetic instruments in order to investigate genetically predicted effects.

We systematically applied the principles of MR to generate effect estimates for each of the 692 proteins on each of the seven neurological phenotypes in turn. This identified 43 genetically predicted effects based on a multiple testing threshold of p < 7.23x10−05 (i.e., 0.05/692 proteins) (Table S3). To support evidence of an effect between proteins and outcomes, we applied genetic colocalization methods to discern whether the causal pQTL at these loci was also responsible for variation in neurological phenotypes.

This identified 12 loci which provided evidence of genetic colocalization based on a PPA > 0.8 for the coloc method or a CLPP > 0.01 for eCAVIAR results (Table S4). Each of these effects has been highlighted on the Manhattan plot in Figure 1. A flowchart illustrating the overall analysis pipeline applied in this study can be found in Figure S1.

Figure 1

A Manhattan Plot to Highlight Genetically Predicted Effects Based on Mendelian Randomization and Genetic Colocalization Analyses on Neurological Phenotypes

Points correspond to the −log10 p values that reflect genetically predicted effects between protein quantitative trait loci and neurological phenotypes. The red dashed line indicates the multiple testing correction applied in analyses (p = 0.05/692 = 7.23 × 10−05). Effects that surpassed this threshold were only included in this plot if they also provided evidence of genetic colocalization, and these effects are colored based on their associated traits.

Genetic Colocalization Helps Develop Insight at Genome-wide Association Loci and Highlights Potentially Novel Signals

Amongst the results with evidence of genetic colocalization were GWAS loci which harbor genes thought to be involved in neurological functionality. For example, the DCC locus, which was associated with neuroticism risk in our analysis (p = 4.29 × 10−11), encodes the netrin-1 receptor, which has been reported to play a role in the development of the nervous system. Similarly, SARM1 was associated with amyotrophic lateral sclerosis risk (p = 1.76 × 10−08) and has been previously reported to influence axonal degeneration. Integrating brain-derived protein data at these GWAS loci can therefore provide insight into the biological pathway by which the underlying causal variants influence these neurological phenotypes.

In contrast, other loci highlighted by our findings are putatively novel trait-associated variants that have not yet reached genome-wide evidence thresholds (i.e., p < 5 × 10−08). As such, evidence of genetic colocalization with pQTL derived from a tissue type relevant to these phenotypes can be valuable in terms of prioritizing loci that are yet to be uncovered by GWAS. Two examples of this include FLOT2 (p = 3.97 × 10−05 with intelligence), which encodes the neuronal signaling factor flotillin-2, and SIDT1 (p = 1.34 × 10−05 with insomnia), which is a dsRNA transporter. Figure 2 illustrates evidence of colocalization at both of these loci with their respective neurological phenotypes and dorsolateral prefrontal cortex proteins. Whereas these signals are likely to be uncovered by GWAS once sample sizes increase, evidence that they colocalize with disease-relevant protein data may shed light on the causative pathway responsible for these signals. We postulate that future endeavors adopting a similar approach but using larger omic datasets will possess increased power to elucidate putatively novel findings.

Figure 2

Locuszoom Plots to Illustrate Evidence of Genetic Colocalization between Proteins and Neurological Phenotypes

Regional −log10 p values at the FLOT2 locus on (A) flottilin-2 levels and (B) intelligence and also at the SIDT1 locus on (C) SIDT1 protein levels and (D) insomnia.

Prioritizing Therapeutic Targets by Undertaking a Phenome-wide Mendelian Randomization Study

For each of the 12 loci highlighted by our initial analysis, we conducted a phenome-wide MR analysis to evaluate putative pleiotropic effects on 700 complex traits and diseases (Table S5). Based on multiple testing comparisons (i.e., p < 0.05/700 = 7.14 × 10−05), along with the same genetic colocalization thresholds used previously, there was evidence of pleiotropy at various loci. For example, based on these criteria, the lead pQTL at SNX32 colocalized with 12 different phenotypes along with Alzheimer disease in the initial analysis (with p = 1.68 × 10−05). As depicted in Figure 3A, therapeutically targeting this gene to reduce risk of Alzheimer disease is genetically predicted to influence other outcomes in the opposite direction, such as HDL cholesterol levels (p = 3.14 × 10−05) and body fat percentage (p = 2.51 × 10−05). Further evaluations of this target are necessary to discern whether increased genetic liability toward Alzheimer disease risk is responsible for these predicted effects (e.g., lower adiposity). Alternatively, genetic variation at this locus may influence these outcomes separately via alternate biological pathways (also known as “horizontal pleiotropy”), which may make SNX32 less attractive as a therapeutic target.

Figure 3

Phenome-wide Association Plots for (A) SNX32 and (B) SARM1 to Investigate Pleiotropic Effects

Each point on these plots corresponds to the −log10 p values derived using the Wald ratio, which are clustered and colored based on the subcategory of each trait and oriented to reflect the direction of effect with each respective protein. Red dashed lines correspond to the multiple testing correction threshold of p < 0.5/700 = 7.14 × 10−05.

In contrast, there were several genes which did not provide evidence of pleiotropy based on this analysis. For instance, the lead pQTL for SARM1, which had an effect on amyotrophic lateral sclerosis risk in the initial analysis (p = 1.76 × 10−08), did not provide evidence of an effect with any of the 700 outcomes assessed based on multiple testing corrections (Figure 3B). We further explored evidence of potential side effects for SARM1 by comparing the distribution of its MR estimates on all 700 traits with those of 500 randomly selected pQTL. The strongest evidence for a secondary effect potentially overlooked by Bonferroni corrections was on coronary artery disease (Ppermutation = 0.004). However, this effect was not supported by evidence of colocalization (PPA = 12.7% and CLPP = 2.92 × 10−04), and this effect had an FDR of 0.06.

Additionally, the only effect surviving Bonferroni corrections for CTSH, which encodes the cathepsin H protein, were on standing and sitting height (p = 4.28 × 10−05 and p = 1.49 × 10−05 respectively). However, these effects were not supported by evidence of genetic colocalization. Examples such as SARM1 and CTSH should therefore be prioritized as worthwhile candidates for therapeutic intervention, given that a lack of pleiotropic effects from this analysis using human genetics supports their safety and efficacy.

Discussion

We have conducted a study to characterize genetic variants associated with neurological phenotypes by harnessing brain-derived protein data. Under the principles of MR, we identified 43 genetically predicted effects across the genome, and this suggested that there may be a shared genetic architecture between neurological phenotypes and the subset of proteins studied. Applying two different genetic colocalization techniques provided evidence that effects at 12 of these loci were driven by a common causal variant. We next undertook a phenome-wide association study for each of these 12 proteins by applying this approach systematically to 700 complex traits and disease endpoints. Doing so elucidated pleiotropic proteins associated with various outcomes, along with protein targets associated more specifically with their corresponding neurological phenotype as identified in our initial analysis.

The influx of high-dimensional datasets concerning intermediate phenotypes provides an exceptional opportunity to unravel the biological mechanisms responsible for GWAS signals. The tissue type used to capture these molecular signatures has been shown to play an important role in such endeavors. For instance, previous comparisons of quantitative trait loci associated with the same gene target identified a correlation of r2 = 0.70 between brain and whole blood. While this suggests that blood may act as a valid proxy for brain tissue the majority of the time, there may be effects that would potentially be overlooked by not using the most pertinent tissue type for the investigated GWAS trait. As an example of this, the fine-mapped pQTL for ERLIN1 in this study shows no evidence on an effect on this protein in whole blood based on the most comprehensive plasma protein QTL analysis to date by Sun et al. (rs11190393, p = 0.92). In fact, the only fine-mapped pQTL we were able to replicate using this data was for CTSH (rs34593439), although this protein has yet to be linked with Alzheimer disease in whole blood analyses based on evidence from the EpiGraphdb platform (URL located in the Web Resources section).

Among the loci highlighted in our study are various genes which have been previously reported to play a role in brain-related activities. This includes DCC, which is associated with neuroticism in our analyses (p = 4.29 × 10−11) and is responsible for expression of the nectin-1 receptor. Nectin-1 has previously been implicated in various neurological and psychiatric disorders, including schizophrenia and depression. Similarly, FLOT2, which was associated with intelligence (p = 3.97 × 10−05), encodes neuronal signaling factor flotillin-2, and it has been linked previously with autism and related disorders. Elsewhere, PSMB4, which encodes a member of the proteasome B-type family, was associated with depression risk (p = 2.36 × 10−05). Proteasomes have been implicated previously in risk of neurodegenerative disorders., These findings may therefore help shed some light on the causal pathway between trait-associated genetic variants at these loci and their respective phenotypes.

Characterizing GWAS signals using tissue-relevant data can also be valuable for translational purposes such as prioritizing therapeutic targets. In particular, SARM1 (p = 1.76 × 10−08 with amyotrophic lateral sclerosis) and CTSH (p = 5.57 × 10−05 with Alzheimer disease) represent the most promising candidates based on our evaluations. This is because our phenome-wide analyses did not detect strong evidence of pleiotropic effects on non-neurological traits, which may foreshadow adverse side-effects from targeting these genes or their mechanism of action using therapeutics. Evidence from the literature has reported that genetic deletion of SARM1 in mice can block pathological axon degeneration., Therapeutically inhibiting SARM1 may therefore be a putatively viable strategy for treating neurodegenerative diseases characterized by axon loss, such as amyotrophic lateral sclerosis. CTSH has also previously been linked with Alzheimer disease, where its expression in the temporal cortices of late-onset Alzheimer patients was shown to be altered. Although our phenome-wide association study included outcomes which are not clinically relevant (e.g., height), the purpose of this “hypothesis-free” analysis was to prioritize potential targets based on overall pleiotropic effects. For example, a target predicted to influence non-clinically relevant endpoints may still be more attractive than one linked only to the target disease being evaluated. It is likely that the majority of therapeutic targets will result in some type of unanticipated side effect, which means that anyone assessing this based on human genetics should primarily be concerned with evaluating whether predicted adverse effects outweigh any potential benefit.

We found, based on evaluations in the GTEx project across 54 tissue types, that various loci highlighted in our analyses were predominantly expressed in brain tissue (Figures S2–S13). However, in our extended analyses based on a previously meta-analyzed sample of 1,194 individuals, only four of the 12 identified proteins provided evidence of genetic colocalization with gene expression (CTSH, KHK, PSMB4, and SNX32, Table S7). There could be various reasons for this lack of agreement, which has been reported previously by the authors of the brainQTL resource, such as technical artifacts in assays. There may also be biological explanations such as canalization, the phenomenon used to describe the robustness of phenotypic characteristics in the presence of abundant genetic variation and environmental conditions., This lack of evidence for the eight proteins which did not colocalize appeared to be due to the fact that their lead pQTL were not also eQTL in the meta-analyzed dataset. However, based on findings from the eQTLGen consortium (n = 31,684), we did find that 10 of the 12 fine-mapped pQTL for these proteins are strongly associated with their corresponding genes’ expression in whole blood. The two exceptions were DCC and RLBP1, which were not analyzed by eQTLGen, and GTEx evaluations suggested that they may not be strongly expressed in whole blood (Figures S4 and S10). As such, a higher proportion of transcriptomic and proteomic signatures may colocalize once sample sizes of brain-derived molecular datasets increase.

This highlights the key limitation of our study, which is the current sample size of accessible proteome-wide data derived from brain tissue (n = 144 from the brainQTL resource). This limited the number of proteins we were able to instrument using pQTL and also meant we were confined to using single-pQTL instruments. Furthermore, it reduced the overall statistical power of the initial pQTL study, which had downstream implications for our colocalization analysis in terms of the number of signals which met conventional thresholds. Future endeavors which continue to uncover the genetic architecture of the human proteome in disease-relevant tissue types will improve our capability to reliably instrument them under the principle of MR. This will also improve the robustness of evidence from genetic colocalization analyses.

We also found that the two genetic colocalization techniques used in this study (the coloc and eCAVIAR approaches) did not always provide corroborating evidence. For example, the loci highlighted in Figure 2 were identified only when using the eCAVIAR approach, despite these plots providing graphical illustrations of genetic colocalization to support findings. A possible explanation for this is the sensitivity of prior distributions selected for the coloc method, which therefore supports our decision to apply multiple methods which make different assumptions about the underlying genetic architecture of a region. However, although genetic colocalization can help support evidence of causality by reducing the likelihood that LD has influenced findings, it cannot rule out horizontal pleiotropy. This is the phenomenon whereby two traits (e.g., a circulating protein and neurological phenotype) are influenced by the same causal variant but by via two independent biological pathways. Functional follow-up work with the proteins highlighted in this work is therefore necessary in order to robustly investigate this.

In conclusion, the findings from this study can help elucidate the causal pathway for genetic variants associated with neurological phenotypes and also prioritize candidate targets for therapeutic intervention. Future studies which leverage increasingly large-scale molecular datasets derived from disease-relevant tissues will continue to develop insights into the mechanisms linking genetic variation to complex traits and disease.

Declaration of Interests

T.R.G. and C.L.R have previously received research funding from Sanofi. T.G.R. and N.K. have both been previously funded by Sanofi. T.R.G receives research funding from GlaxoSmithKline and Biogen. However, none of these factors contributed to the design or analysis of this study.