ARTICLE https://doi.org/10.1038/s41467-020-14483-x OPEN Regulatory sites for splicing in human basal ganglia are enriched for disease-relevant information Sebastian Guelfi1,72, Karishma D’Sa1,2,72, Juan A. Botía1,3,72, Jana Vandrovcova1, Regina H. Reynolds 1, David Zhang1, Daniah Trabzuni 1,4, Leonardo Collado-Torres 5, Andrew Thomason 6, Pedro Quijada Leyton6, Sarah A. Gagliano Taliun 7, Mike A. Nalls8,9, International Parkinson’s Disease Genomics Consortium (IPDGC)*, UK Brain Expression Consortium (UKBEC)*, Kerrin S. Small 10, Colin Smith11, Adaikalavan Ramasamy1,2,12, John Hardy1, Michael E. Weale2,13,72 & Mina Ryten1,72 ✉ 1234567890():,; Genome-wide association studies have generated an increasing number of common genetic variants associated with neurological and psychiatric disease risk. An improved under- standing of the genetic control of gene expression in human brain is vital considering this is the likely modus operandum for many causal variants. However, human brain sampling complexities limit the explanatory power of brain-related expression quantitative trait loci (eQTL) and allele-specific expression (ASE) signals. We address this, using paired genomic and transcriptomic data from putamen and substantia nigra from 117 human brains, inter- rogating regulation at different RNA processing stages and uncovering novel transcripts. We identify disease-relevant regulatory loci, find that splicing eQTLs are enriched for regulatory information of neuron-specific genes, that ASEs provide cell-specific regulatory information with evidence for cellular specificity, and that incomplete annotation of the brain tran- scriptome limits interpretation of risk loci for neuropsychiatric disease. This resource of regulatory data is accessible through our web server, http://braineacv2.inf.um.es/. 1 Reta Lila Weston Research Laboratories, Department of Molecular Neuroscience, University College London (UCL) Institute of Neurology, London, UK. 2 Department of Medical & Molecular Genetics, School of Medical Sciences, King’s College London, Guy’s Hospital, London, UK. 3 Departamento de Ingeniería de la Información y las Comunicaciones, Universidad de Murcia, Murcia, Spain. 4 Department of Genetics, King Faisal Specialist Hospital and Research Centre, Riyadh, Saudi Arabia. 5 Lieber Institute for Brain Development, 855 North Wolfe Street, Baltimore, MD, USA. 6 Goldsmiths, University of London, New Cross, London, UK. 7 Center for Statistical Genetics, University of Michigan, Ann Arbor, MI, USA. 8 Laboratory of Neurogenetics, National Institute on Aging, US National Institutes of Health, Bethesda, MD, USA. 9 Data Tecnica International, Glen Echo, MD, USA. 10 Department of Twin Research and Genetic Epidemiology, King’s College London, London, UK. 11 Department of Neuropathology, MRC Sudden Death Brain Bank Project, University of Edinburgh, Edinburgh, UK. 12 Singapore Institute for Clinical Sciences, Brenner Centre for Molecular Medicine, Singapore, Singapore. 13 Genomics plc, Oxford, UK. 72These authors contributed equally: Sebastian Guelfi, Karishma D’Sa, Juan A. Botía, Michael E. Weale, Mina Ryten. *A full list of consortia members appears at the end of the paper. ✉email: mina.ryten@ucl.ac.uk NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x T he use of genome-wide genotyping in large patient and these paired data, we searched for eQTL associations between control populations has resulted in the identification of ~6.5 million genetic variants and ~411,000 RNA expression traits increasing numbers of common variants that impact on in putamen and ~370,000 RNA expression traits in substantia the risk of a wide range of neurological and psychiatric condi- nigra, resulting in ~5.3 billion eQTL tests. We generated RNA tions, including Parkinson’s disease1–3, Alzheimer’s disease4–6, expression traits using both annotation-based (known tran- and schizophrenia7,8. However, the majority of these risk loci are scripts) and annotation-agnostic approaches, noting that both still poorly characterised, and we do not yet fully understand the types of traits distinguished between the brain regions (Fig. 1a, underlying molecular and cellular processes through which they Supplementary Fig. 2). Within annotated regions, RNA quanti- act. As it is reasonable to assume that many causal variants fication was performed with RNA processing in mind (Fig. 1a) to operate by regulating gene expression, several studies have produce four separate measures of transcription, which formed attempted to address this problem through the use of expression the bases of our eQTL analysis. This resulted in the generation of quantitative trait loci (eQTL) and allele-specific expression (ASE) four types of eQTLs (gi-eQTLs, e-eQTLs, ex-ex-eQTLs, and ge- analyses in a wide range of human tissues, with the aim of finding eQTLs), of which two were designed to capture the genetic reg- eQTL and ASE signals that colocalise with disease risk signals9,10. ulation of splicing eQTLs (e-eQTLs and ex-ex-eQTLs). Finally, This approach has had success, but perhaps not as much as we included annotation-independent approaches to quantify might have been expected for all diseases. Although the identi- transcription. We focused specifically on unannotated transcrip- fication of eQTLs in blood has provided insights into auto- tion within intergenic regions (producing i-eQTLs, Methods). immune disorders11,12, the utility of brain-related eQTL and ASE Following stepwise conditional analyses under a false discovery data sets, particularly with regard to neurodegenerative disorders, rate (FDR) of 5%, we identified 19,156 separate significant eQTL has been harder to demonstrate. For example, monocyte eQTL signals (hereafter, eQTLs) genome-wide (Supplementary Data 2–6), data sets appear to provide greater insights for Alzheimer’s dis- of which 359 were secondary eQTLs (i.e., eQTLs with independent ease13, probably because they reflect regulatory processes in effects after conditioning on the primary eQTL in the region). microglia. This would suggest that while relevant eQTL and ASE Although there was a substantial difference in the number of signals are present in the brain, they are currently difficult to eQTLs identified in putamen and substantia nigra, the difference in detect given the constraints on eQTL and ASE analyses in sample size (N = 111 for putamen and N = 69 for substantia nigra) human brain. most likely accounted for this. However, eQTL discovery was not At present, the most easily available sampling method for brain simply driven by the number of features tested. Notably, the rate of tissue is post mortem, making repeat sampling impossible and eQTL discovery, defined as the percentage of all expression features typically leading to smaller sample sizes, particularly for smaller tested with at least one significant associated eQTL, was highest in structures such as the substantia nigra. Furthermore, the brain is a unannotated intronic and intergenic regions (Fig. 1b), suggesting highly complex organ. Not only is it split into many regions with that such eQTLs could be biologically important. known inter-regional differences in expression14, each region is composed of an assemblage of different cell types, which com- plicates the interpretation of transcriptomic data and limits sta- eQTLs and i-eQTL target regions show high replication rates. tistical power. Finally, the brain transcriptome is unusual in We found that 50.6 and 50.4% of the testable eQTLs identified in having a high degree of alternative splicing and a high degree of putamen and substantia nigra, respectively, were also detected non-coding RNA activity15, much of which has yet to be fully using microarray data, generated by the UK Brain Expression characterised16. Consortium17 and based on a common set of RNA samples. We We address the latter of these constraints by conducting total also found that 39.3% of putamen and 50.6% of substantia nigra RNA sequencing in two basal ganglia regions of clinical interest to eQTLs replicated in the GTEx (v7) data resource19, using their human neurodegenerative and neuropsychiatric disorders: the 111 putamen and 80 substantia nigra samples respectively. Fur- substantia nigra and putamen. Using a comprehensive set of thermore, we investigated eQTL replication across all brain analyses to interrogate different stages of RNA processing and regions studied in GTEx (Supplementary Table 1). We demon- uncover novel unannotated transcripts, we seek to identify not strated a replication rate of 19 to 53.3% for putamen, with the only disease-relevant regulatory loci but also the types of analyses highest replication rates observed in cerebellum (53.3%) and and regulatory positions yielding the most brain and disease- caudate (48.2%). In the case of substantia nigra samples, repli- specific information. We find that splicing eQTLs are enriched for cation rates were more similar across all brain tissues neuron-specific regulatory information; that ASE analyses, prob- (50.6–62.0%), potentially reflecting the relatively low sample ably by more effectively controlling for cellular heterogeneity, numbers used in the eQTL analysis in both our study (N = 65) provide highly cell-specific regulatory information; and that and that performed by GTEx (N = 80). Furthermore, we inves- incomplete annotation of the brain transcriptome is limiting the tigated replication using eQTLs from the larger dorsolateral interpretation of risk loci for neuropsychiatric disease. We release prefrontal cortex data sets generated by the PsychENCODE and the rich resource of eQTL and ASE data generated in this study CommonMind consortia20,21. We found that although the eQTL through a searchable web server, http://braineacv2.inf.um.es/ data sets generated by PsychENCODE and CommonMind con- (Supplementary Fig. 1). sortia were based on the analysis of sample sets, which are ~13 and ~5 times greater in size than GTEx, the replication rates were (64.2% and 54.2%, respectively, as compared with 52.7% when Results using the basal ganglia GTEx tissues). In contrast, when we RNA quantification and eQTL discovery. We assayed DNA and checked our eQTL signals against those reported by Lappalainen RNA from 180 brain samples originating from 117 individuals of and colleagues10 using RNA-seq analysis of 373 lymphoblastoid European descent, which were part of the UK Brain Expression cell lines, we found that despite the larger sample size in this Consortium data set17 and which were classified as neurologically study, only 22.0% of putamen and 24.2% of substantia nigra healthy based on the absence of neurological disease during life eQTLs were replicated. and neuropathological assessment (Supplementary Data 1). We Given the paucity of existing eQTL analyses using annotation- focused on putamen and substantia nigra samples owing to their independent approaches, we focused on validating the expression distinctive expression profiles and disease relevance18,19. Using of unannotated intergenic regions that were the target of a 2 NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x ARTICLE a RNA type Pre-mRNA Transcript-specific Transcribed intergenic Genic Gene Transcript Gene A Gene B Gene Read alignment Read counts Gene-intronic counts Exon-exon Exon counts Intergenic counts Gene-exonic counts junction counts Quantification In-house Altrans DEXseq Derfinder Htseq-count method pipeline eQTL type gi-eQTLs ex-ex-eQTLs e-eQTLs i-eQTLs ge-eQTLs Biological Transcriptional rate Alternative splicing Unannotated expression Steady-state mRNA relevance b 6 5 4 Yield (%) 3 2 1 0 gi-eQTLs ex-ex-eQTLs e-eQTLs i-eQTLs ge-eQTLs eQTL type Fig. 1 Similar eQTL yield for unannotated expression features compared with annotated features. a Overview of transcriptome quantification. RNA was quantified using five pipelines, each targeting distinct stages of RNA processing, and each followed by eQTL generation. Within annotated regions of the transcriptome, reads were mapped to expression features and thereafter RNA was quantified. These features included: all intronic and exonic regions of a gene (producing gene-intronic gi-eQTLs and gene-exonic ge-eQTLs, respectively); individual exons (producing e-eQTLs); and exon–exon junctions (producing ex-ex-eQTLs). As total RNA was used for library construction, reads mapping to introns were presumed to be owing to pre-mRNA within samples (an assumption supported by previous analyses using a subset of these data69). Quantification of individual exons and exon–exon junctions provided a means of identifying loci that impact on alternative splicing. In common with most eQTL analyses, we also calculated overall gene expression using all reads mapping to exons of a given gene, resulting in an expression metric that is influenced by transcriptional rate, splicing and RNA degradation rates. Finally, we included annotation-independent approaches to quantify transcription. We focused specifically on unannotated transcription within intergenic regions (producing i-eQTLs, Online Methods). b eQTL yields for both tissues were calculated as the number of expression features within a category with at least one significantly associated eQTL divided by the total number of tested features within the same category. Source data are provided as a Source Data file. significant i-eQTL. Using data provided by GTEx and processed in overlap with enhancer regions as defined within the for re-use by recount222, we found that 70.3% of all such i-eQTL GeneHancer database v4.423 among i-eQTL target regions target regions (in putamen and substantia nigra combined) were (17.9%) as compared with e-eQTL targets (5.9%), suggesting that detected in at least one other human tissue within the GTEx data the transcribed regions targeted by i-eQTLs are highly enriched set, with the highest validation rates observed among brain for eRNAs. regions (Supplementary Table 1). We also explored the possibility We further characterised i-eQTL target regions based on their that the transcribed regions detected in our analysis and regulated relationship to known genes (Fig. 2a). Using reads spanning by i-eQTLs may represent enhancer RNAs as another means of known exons and novel regions, physical proximity and understanding the biological relevance of our findings. Consider- correlation in expression, we categorised unannotated expressed ing all transcribed regions targeted by an i-eQTL and accounting regions into those with strong, moderate, or weak evidence for for their genomic size, we found that there was a 3.0-fold increase being part of a known gene (Fig. 2a, Online Methods). This NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x a b i-eQTL target region Strong Moderate Weak Uncharacterised Categorisation features Output categories Biological interpretation 1.0 i-eQTL target 0.8 R2 with reference gene 0.6 Exon junction? Yes Strong 0.4 No 0.2 High co- expression and < 5 kb Yes Moderate 0.0 to nearest gene? 0 10 100 1K 5K 100 K 1M Distance from reference gene (bp) No c 100 Low co- PUTM & SNIG expression and 75 > 5 kb from nearest Yes Weak Validation (%) gene? Brain 50 Strong No Moderate Weak 25 Uncharacterised Uncharacterised 0 Strong Moderate Weak Uncharacterised i-eQTL target region d Strong Moderate Weak Chr 10 Chr 7 Chr 6 Chromosome 13.179 mb 13.181 mb 13.183 mb 13.185 mb 13.187 mb 50.515 mb 50.516 mb 50.517 mb 50.518 mb 50.519 mb 28.389 mb 28.390 mb 28.391 mb 28.392 mb 28.393 mb Primers Split read Sanger sequence Intergenic region Gene track OPTN FIGNL1 Fig. 2 i-eQTL target regions show high replication in independent data sets and validate experimentally. a Characterisation of i-eQTL target regions (unannotated expressed regions that were the target of a significant i-eQTL) was based on several features reflecting their relationship to known genes. These features were used to classify these regions into those with strong, moderate, and weak evidence for being part of a known gene. Regions categorised as strong and moderate are considered likely to be novel exons of known genes or misannotations of existing exon boundaries, whereas weak regions are presumed to be independent of any known genes. b Scatterplot of genomic distance and correlation of expression between i-eQTL target regions and their reference genes. c The expression of unannotated expressed regions was validated in GTEx data, using brain region-specific and global brain expression data. Validation rates in putamen and substantia nigra GTEx expression data were combined and displayed separately from validation rates in RNA-seq data from all GTEx brain regions. d Sequencing results for i-eQTL target regions with strong, moderate, and weak evidence of being part of a gene. In each case, tracks are provided relating to the location of the primers used to amplify the unannotated expressed region, the RNA-seq split read, the alignment of Sanger-sequenced cDNA, and the predicted boundaries of the unannotated expression region. Source data are provided as a Source Data file. approach allowed us to characterise 68.1% of all unannotated i-eQTLs and non-standard eQTLs are largely novel signals. As expression regions (Fig. 2b). The validation rate for expression in 51.2% of our characterised i-eQTL target regions have strong or the GTEx data resource was 98.5% for unannotated expressed moderate evidence linking them to a known gene, we further regions with strong evidence for being part of a known gene, classified i-eQTLs into those with evidence for being new reg- 93.8% for those with moderate evidence, and still high at 60.4% ulatory variants versus those appearing to act in a consistent for regions with weak evidence for being part of a known gene manner across all exons (thus recapitulating gene-level signals). (Fig. 2c). We also selected eight unannotated expressed regions Using a modified test of heterogeneity (Methods) we separately for experimental validation (Supplementary Table 2). Regardless analysed i-eQTLs with strong, moderate, and weak evidence for of their putative relationship to existing genes, all eight regions being linked to a known gene (Fig. 3). This analysis demonstrated validated using Sanger sequencing (Fig. 2d). In the case of that some i-eQTLs were indeed re-discovered versions of existing unannotated expressed regions with moderate evidence for eQTL signals (Fig. 3a). However, many i-eQTLs appear to be association, this analysis also enabled us to clarify the exon independent regulatory sites. For example, SNP rs4696709 reg- boundaries. For example, sequencing confirmed the existence of a ulates DER10633 expression, a probable novel exon of ABLIM2 novel exon of FIGNL1 (DER18381, Fig. 2d). Thus, using a (based on the presence of junction reads), but there is no sig- combination of public data resources and experimental work, we nificant co-regulation of other exons of ABLIM2 (Fig. 3b). Even demonstrated the validity of annotation-independent approaches among those i-eQTLs with strong evidence linking them to a in transcriptomic analysis. known gene, the percentage of i-eQTLs sharing signals with 4 NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x ARTICLE a Chromosome 13 43.6 mb 43.62 mb 43.64 mb 43.66 mb 43.68 mb 4 43.61 mb 43.63 mb 43.65 mb 43.67 mb 43.69 mb PEER corrected expression 1 –log10 (FDR) i-eQTL 3 2 0 1 –1 0 4 GG TG TT Genotype –log10 (FDR) 3 ge-eQTL 2 PEER corrected expression 0.5 1 0 0.0 –0.5 Gene track DNAJC15 –1.0 GG TG TT i-eQTL region Genotype b Chromosome 4 7.95 mb 7.97 mb 7.99 mb 8.01 mb 7.96 mb 7.98 mb 8 mb 4 PEER corrected expression 1 –log10 (FDR) 3 i-eQTL 2 0 1 –1 0 4 AA GA GG Genotype –log10 (FDR) 3 ge-eQTL 2 PEER corrected expression 0.5 1 0.0 0 i-eQTL region –0.5 Gene track –1.0 ABLIM2 AA GA GG Genotype c 100 d 100 Distinct eQTL Shared eQTL Distinct eQTL Shared eQTL 75 75 % of eQTL signals % of regions 50 50 gi-eQTL Strong ex-ex-eQTL Moderate e-eQTL Weak 25 25 i-eQTL 0 0 Strong Moderate Weak gi-eQTL ex-ex-eQTL e-eQTL i-eQTL eQTL type eQTL type known annotation expression features was only 44% (Fig. 3c). We that the regulatory information captured by these eQTLs is extended this analysis to determine whether i-eQTLs with strong indeed novel. Thus, across all types of characterised i-eQTLs, we or moderate evidence linking them to a known gene could found evidence for the majority representing novel regulatory nonetheless be found in the larger eQTL data sets provided by variants, acting in a transcript-specific manner. CommonMind, PsychENCODE, and BrainSeq20,21,24. In fact, the We also asked whether our alternative annotation-based eQTL majority of i-eQTLs were not identified (CommonMind, 73.0%; classes (gi-eQTLs, e-eQTLs, and ex-ex-eQTLs) provided novel PsychENCODE, 82.1%, and BrainSeq, 70.0%) again suggesting regulatory information compared with the standard gene-level NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x Fig. 3 i-eQTL target regions have evidence for distinct regulation. a Local association plots (−log10 FDR-corrected p values for eQTL association), illustrating sharing of the rs113317084 variant (red point) between the i-eQTL-targeted region, DER32583 (green track), and the ge-eQTL-targeted gene, DNAJC15 (blue track). b Local association plot illustrating no sharing of the rs4696709 variant (red point) between the i-eQTL-targeted region, DER10633 (green track), and the ge-eQTL-targeted gene, ABLIM2 (blue track). The detection of reads spanning DER10633 and an annotated exon within ABLIM2 provides compelling evidence that this region represents a novel exon of the gene. c Heterogeneity (distinct vs. shared) of i-eQTL signals, cross- categorised by the strength of evidence linking their target region to a known gene, suggests that most are distinct and likely represent novel regulatory variants acting in a transcript-specific manner. Heterogeneity was determined using a modified beta-heterogeneity test, accounting for the dependency structure arising from within-individual and within-gene correlations. i-eQTL beta-coefficients were compared with that of the known exon with most evidence of association with the i-eQTL target region. All eQTL signals with an FDR-corrected p value for heterogeneity < 0.05 were considered distinct, whereas those with an FDR-corrected p value > 0.05 were considered shared (similar beta-coefficients). d Heterogeneity (distinct vs. shared) of non- standard eQTL classes (gi-eQTLs, e-eQTLs, ex-ex-eQTLs, and i-eQTLs) suggests that many of these classes are distinctly regulated. Heterogeneity was determined using a modified beta-heterogeneity test comparing beta-coefficients from ge-eQTLs to those derived from non-standard eQTL analyses applied to the same gene. This analysis was performed separately for gi-eQTLs (tagging pre-mRNA), e-eQTLs, and ex-ex-eQTLs (tagging splicing) and all i- eQTLs (tagging unannotated expression). All eQTL signals with an FDR-corrected p value < 0.05 were considered distinct, whereas an FDR-corrected p value > 0.05 was taken as evidence of eQTL sharing. Source data are provided as a Source Data file. eQTL analysis (ge-eQTLs). Again we used a modified test of beta- adjusting for the size of each annotation, we found that unan- heterogeneity to determine eQTL signal sharing among these notated transcribed regions and exons had a similar level of classes. Although 66.8% of gi-eQTLs were detectable with overlap with brain-relevant risk loci (3.5% for exons and 2.9% for standard ge-eQTLs, this figure was only 6.9% for splicing eQTLs transcribed intergenic regions after adjustment for annotation (Fig. 3d, suggesting that our additional eQTL classes provided size). However, the enrichment of brain-relevant risk loci was distinct regulatory information driven by splicing effects. higher for novel transcribed regions as compared to exons (1.51- fold for targets of e-eQTLs versus 2.70 fold for targets of i-eQTLs Splicing eQTLs are enriched for neuronal information. We after adjustment for annotation size). Furthermore, we found a asked whether our different eQTL classes varied in terms of the significant enrichment for GWAS variants that were associated cellular specificity of their target expression features. We used with neurological and behavioural disorders compared with all weighted gene co-expression network analysis, in combination with other SNP-phenotype associations (Supplementary Fig. 3) among publicly available cell-specific annotation data, to assign eQTL our eQTLs. Although the enrichment of brain-relevant GWAS target expression features to one of five broad cell types: neuron, associations was most evident in ex-ex-eQTLs and gi-eQTLs, we oligodendrocyte, astrocyte, microglia, and endothelial cell (Fig. 4a, also found a significant enrichment for i-eQTLs (FDR-corrected Online Methods). This module membership approach allowed us to Fisher’s Exact p value = 6.45 × 10−7). i-eQTLs provided useful provide putative cellular classifications for expression features even information for 36.7% of all the neurologically relevant risk loci if they were outside known annotations. We confidently assigned within this analysis (equating to 76 loci). Given that these find- up to 75% of all analysed genes to a specific cell type, and these were ings could potentially be driven by correlations between i-eQTLs then related to 41.5% of all eQTL target expression features. We and more conventional eQTL signals, we repeated this analysis observed a significant enrichment of neuronal genes in all non- only using i-eQTLs considered to be independent regulatory standard eQTL classes in one or both tissues investigated (Fig. 4b, signals (based on modified beta-heterogeneity testing described Supplementary Table 2). These included i-eQTLs targeting unan- above). We found that the enrichment of brain-relevant risk loci notated expressed regions (FDR-corrected Fisher’s Exact p value = among i-eQTLs increased in significance in this sub-group (FDR- 1.20 × 10−2 in putamen, Fig. 4b). Furthermore, we found that the corrected Fisher’s Exact p value = 7.01 × 10−8). Thus, our analysis targets of splicing eQTLs were significantly enriched for neuronal suggests that i-eQTLs do contribute to the understanding of a genes (FDR-corrected Fisher’s Exact p values = 1.21 × 10−7 and significant proportion of neurologically relevant risk loci. 2.28 × 10−5 for e-eQTLs and ex-ex-eQTLs, respectively, in sub- We further explored signal enrichment in i-eQTLs in relation to stantia nigra), oligodendrocyte genes (FDR-corrected Fisher’s Exact two neurological diseases related to basal ganglia dysfunction: p value = 1.7 × 10−3 and 4.1 × 10−2 for e-eQTLs and ex-ex-eQTLs, Parkinson’s disease and schizophrenia. Using GWAS summary respectively in substantia nigra) and astrocyte genes (FDR-corrected statistics for these diseases1,7, we performed colocalisation analyses Fisher’s Exact p value = 8.74 × 10−4 and 1.12 × 10−3 for e-eQTLs for disease-risk association signals against i-eQTL signals using the and ex-ex-eQTLs respectively in substantia nigra, Fig. 4b). This coloc26 programme. We identified 23 i-eQTL signals that colocalised points to the importance of capturing splicing information in the with risk loci for schizophrenia or Parkinson’s disease (Supple- analysis of human brain samples. mentary Data 8). Among the former, we identified a signal indexed by the lead SNP rs35774874 (GWAS p value = 1.97 × 10−11) that i-eQTLs are enriched for disease-relevant information. We colocalised with an i-eQTL targeting a probable novel 3′-UTR of investigated the overlap of unannotated transcribed regions and SNX19 (posterior colocalisation probability = 0.75), a gene that has eQTL sites with known GWAS association signals. We used the already been highlighted in schizophrenia27,28. Similarly, we US National Human Genome Research Institute/European identified a colocalisation of the Parkinson’s disease GWAS lead Bioinformatics Institute (NHGRI-EBI) GWAS catalogue, restric- SNP rs4566208 (GWAS p value = 2.28 × 10−7) with an i-eQTL ted to genome-wide significant SNPs (P < 5 × 10−8) and stratify- regulating a probable novel exon of ZSWIM7 (i-eQTL p value = ing for SNP-phenotype associations of relevance to neurological/ 1.09 × 10−5; posterior colocalisation probability = 0.89, Supplemen- behavioural disorders as defined within the STOPGAP data- tary Fig. 4). However, we also found seven co-localising i-eQTL base25. First, we investigated the possibility that unannotated signals targeting unannotated expressed regions that were not transcribed regions could themselves harbour risk loci of rele- linked to a known gene. For example, we found that the vance to neurological diseases, by calculating the proportion of schizophrenia risk SNP rs12908161 (GWAS p value = 9.41 × transcribed intergenic regions containing a brain-relevant risk 10−10) had a posterior colocalisation probability of 1.00 with an locus and comparing this value to that for expressed exons. After i-eQTL targeting the unannotated expressed region DER36302 6 NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x ARTICLE a Cell type-specific markers Gene expression Co-expression networks Assign modules to WGCNA cell types based on Cell type-specific with k-means enrichment Modules Filter for genes with MM > 0.3 eQTLs Filter for genes/ERs Fisher’s exact test Cell type-specific eQTL target Cell type-specific with MM > 0.3 for enrichment enrichments genes/ERs eQTL genes/ERs b Substantia nigra Putamen gi-eQTLs gi-eQTLs Transcriptional rate ex-ex-eQTLs ex-ex-eQTLs Alternative splicing Cell type e-eQTLs e-eQTLs Alternative Astrocyte splicing Endothelial Microglia Neuron Oligodendrocyte i-eQTLs i-eQTLs Unannotated expression ge-eQTLs ge-eQTLs Steady-state mRNA 0 2 4 6 0 5 10 15 20 25 –log10 (p value) –log10 (p value) Fig. 4 Non-standard eQTL analyses produce additional biologically relevant information. a Schematic diagram showing the use of gene co-expression networks to assign eQTL target genes and unannotated expressed regions (ERs) to the cell type most likely to be driving gene expression in the tissue. We used the WGCNA R package63. b eQTL classes were variably enriched for genes with cell-biased expression, highlighting the importance of capturing this information. Enrichment of genes with cell-biased expression within eQTL targeted expression features was performed separately for each tissue and was determined using a Fisher’s Exact test and a significance cutoff of P < 0.05 (dashed red line at −log10(P) = 1.30). Genes assigned to modules significantly enriched for brain-related cell type markers and with a module membership of > 0.3 were allocated a cell type. Next, for each eQTL targeting a known genic region or an unannotated expressed region with high or moderate evidence linking it to a known gene, if the target gene was allocated to a cell type then the related eQTL received the same cell type label. For eQTLs targeting unannotated expressed regions with low evidence for association with a known gene or which could not be classified, we assigned the target expression feature to a module (and by inference a cell type) based on its highest module membership providing the module membership was at least 0.3. Finally, for each eQTL class and each cell type, namely neuron, microglia, astrocyte, oligodendrocyte, and endothelial cell, we applied a Fisher’s Exact test to test for enrichment of that cell type label among the genes associated to the eQTL class. Expression features targeted by different eQTL classes were variably enriched for genes with cell-biased expression, highlighting the importance of capturing this information. Enrichment of genes with cell-biased expression within eQTL targeted expression features was performed separately for each tissue and was determined using a Fisher’s Exact test and a significance cutoff of P < 0.05 (dashed red line at −log10(P) = 1.30). Source data are provided as a Source Data file. (chr15:84833811-84833975, eQTL p value = 7.94 × 10−10 in puta- ASE discovery and validation. We applied ASE analysis to a men, Fig. 5a). This novel transcribed region appears to be subset of 84 brain samples (substantia nigra n = 35; putamen n = independent of neighbouring genes, and is expressed in human 49) for which we had access to whole-exome sequencing in brain, with the highest expression in the anterior cingulate and addition to SNP genotyping data (Methods, Supplementary frontal cortex, two brain regions relevant to schizophrenia (Fig. 5b). Data 1). ASE analysis quantifies the variation in expression NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x a Type SCZ GWAS eQTL b Tissues Brain Non brain 80 60 RPKM 7.5 40 –log10 (p value) 20 5.0 0 Brain-Anterior cingulate cortex (BA24) Brain-Frontal cortex(BA9) Brain-Hypothalamus Brain-Nucleus accumbens (basal ganglia) Brain-Hippocampus Brain-caudate (basal ganglia) Brain-Putamen (basal ganglia) Brain-Amygdala Pituitary Brain-Substantia nigra Brain-Spinal cord (Cervical c-1) Bladder Colon-Sigmoid Nerve-Tibial Esophagus-Muscularis Esophagus-Gastroesophageal junction Heart-Atrial appendage Adipose-Subcutaneous Skin-Sun exposed (Lower leg) Brain-Cerebellar hemisphere Colon-Transverse Artery-Coronary Lung Skin-Not sun exposed (Suprapubic) Artery-Aorta Small intestine-Terminal lleum Adrenal gland Artery-Tibial Adipose-Visceral (Omentum) Spleen Stomach Minor salivary gland Esophagus-Mucosa Thyroid Heart-Left ventricle Pancreas Kidney-Cortex Muscle-Skeletal Liver Whole blood 2.5 0.0 84.0 84.5 85.0 85.5 Tissue Position (Mb) Fig. 5 Annotation-independent approaches yield disease-relevant information. a Colocalisation of the schizophrenia GWAS lead SNP rs950169 (GWAS p value = 7.62 × 10−11) and the i-eQTL targeting DER36302 (eQTL p value = 1.15 × 10−10 in putamen). b Expression of DER36302 across tissues sampled by the GTEx consortium. Brain tissues are highlighted in yellow with the anterior cingulate and frontal cortex showing the highest expression. Source data are provided as a Source Data file. a Imprinting Nonsense-mediated decay cis-eQTL T T G T C C A C T T T T T T T T T T C T T C C C C C b 40 c LMBRD2 PEER corrected expression 30 Consistent 0.2 % ASE signals 20 0.0 Inconsistent 10 –0.2 0 0.0 0.2 0.4 0.6 0.8 1.0 CC CA Proportion of positive ASE signals Genotype amongst individuals Fig. 6 Allele-specific expression provides evidence of dosage compensation. a Overview of the mechanisms by which allele-specific expression can arise. Allele-specific expression can arise through epigenetic effects (e.g., imprinting), heterozygous mutations triggering nonsense-mediated decay of transcripts, and regulation by (for example) a cis-regulatory variant (cis-eQTL). b The majority of allele-specific expression signals passing FDR < 0.05 produced unidirectional signals (0 or 1) and were considered consistent. Inconsistent ASE signals (those that were not unidirectional in ≥ 10 individuals) were found only in known imprinted genes, thus providing additional validation of our ASE signals. c Comparison of LMBRD2 expression in putamen from one individual heterozygous for a rare stop gain mutation (CA) in the gene versus all other individuals (CC) revealed a significant reduction in LMBRD2 expression, implying effective nonsense-mediated decay. Data presented using Tukey-style box plots. Source data are provided as a Source Data file. between two haplotypes of a diploid individual distinguished by found in putamen and 11,871 in substantia nigra (Supplementary heterozygous genetic variation, and so can capture the effects of a Data 9). Consistent with previous studies, we found that ASEs can range of regulatory processes, namely genomic imprinting, operate as markers of imprinting or parent-of-origin effects: ASE nonsense-mediated decay and cis-regulation (Fig. 6a). In total signals that are not unidirectional across individuals are expected 252,742 valid heterozygous SNPs (hetSNPs) across 53 individuals to be enriched for imprinted genes (Fig. 6a). Consistent with were analysed. Of these, 7.62% (19,266) were significant ASE expectation, of all genes containing an ASE, 170 were identified signals (hereafter, ASEs) at FDR < 5% in at least one sample, on the X chromosome (equating to 1.96%). Furthermore, we covering 8654 genes. Of the 19,266 ASEs identified, 12,096 were observed that all inconsistent ASE signals (those that were not 8 NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x ARTICLE unidirectional within ≥ 10 individuals) were located within genes (Fisher’s Exact test p value = 4.50 × 10−5 and p value = 2.07 × known to be imprinted (as reported in www.geneimprint.com or 10−19 using a randomisation approach). This strongly suggests within the literature27–30) compared with 1–5% of consistent that ASEs are enriched for variants with effects on splicing. signals (Fig. 6b). We also found evidence for the generation of ASE signals ASEs show biologically and disease-relevant enrichments. We through nonsense-mediated decay. We identified 61 protein- assessed the cellular specificity of the regulatory information truncating variants (defined as stop gain, donor splice site, and ASEs provide. Given that ASE analysis is performed within an donor acceptor mutations) among our ASEs. Consistent with individual and so is not subject to the confounding effects of expectation, the majority of these variants were predicted to cause cellular heterogeneity across individuals, we expected that ASEs nonsense-mediated decay (52.5% using SNPeff31) and appeared would be a powerful means of obtaining cell-specific regulatory to result in mono-allelic expression, with >95% of all reads at the information. Using a similar approach to that applied to eQTLs to ASE site originating from a single allele. These extreme ASEs are assign genes containing ASEs to brain-relevant cell types (neu- expected to generate an effective reduction in gene dosage, and rons, oligodendrocytes, astrocytes, microglia, and endothelial therefore to cause a significant reduction in the total expression of cells), we found that ASE-containing genes were highly enriched the affected genes. An example of this pattern is seen in the for neuronally expressed genes (Fisher’s Exact test FDR-corrected LMBRD2 gene (Fig. 6c). p values of 9.97 × 10−235 in putamen and 3.05 × 10−97 in sub- Finally, to check the overall reliability of our findings, we stantia nigra). We also found significant enrichments (FDR looked for validation of our ASEs in an independent data set of p value < 5%) for oligodendrocyte, astrocyte, microglia, and 462 lymphoblastoid cell lines reported by Lappalainen and endothelial gene sets. Although we observed a similar pattern of colleagues10. We found that 67% of testable ASE signals could be cell type-specific enrichments among eQTLs, the strength of detected at an FDR < 5%, demonstrating the reliability of our ASE evidence for cellular specificity of ASEs was striking, suggesting sites while also suggesting the presence of brain-specific ASEs. that incomplete covariate correction may be hampering the power of eQTL analyses (Fig. 7a, Supplementary Data 7). Finally, we used GWAS summary data sets for Parkinson’s ASEs tag both gene level and splicing eQTLs. Cis-regulatory disease1 and schizophrenia7 to investigate the disease relevance of variants are known to be one important generator of ASEs32. We ASEs. As GWAS loci often lie close to genic regions and so are therefore investigated the overlap of eQTLs with ASEs in our likely to overlap by chance with ASE signals, we used a data, and compared it with the overlap observed with randomly randomisation approach34 to investigate the enrichment of selected non-ASE heterozygous SNPs. After controlling for read GWAS loci within our ASEs (Methods). We compared overlaps depth, we identified a highly significant enrichment of eQTLs between risk loci and ASEs to overlaps between risk loci and among our ASEs (p value = 2.65 × 10−195 in putamen and 9.99 × randomly selected non-ASE sites, and found a highly significant 10−111 in substantia nigra, using a randomisation approach—see enrichment of GWAS risk loci for both schizophrenia (p value = Online Methods). This enrichment remained significant when we 7.49 × 10−35, using a randomisation approach, for ASEs derived restricted our analysis to eQTLs with effects on splicing (e-eQTLs from both tissues) and Parkinson’s disease (p value = 4.19 × 10−7, and ex-ex-eQTLs, p value = 1.17 × 10−178 in putamen and 1.29 × using a randomisation approach, for ASEs derived from both 10−89 in substantia nigra, using a randomisation approach). tissues). We validated these findings using stratified LD score However, as expected, it was absent when we considered ASEs regression35 by treating our ASE sites as a form of binary located within imprinted genes, where the parental origin of the annotation, and interestingly found that the enrichment of SNP rather than the impact of cis-regulatory sites is expected to Parkinson’s disease heritability appeared to be more specific to drive allele-specific expression (p value = 0.923 in putamen and ASEs identified in substantia nigra using this approach (Fig. 7b, 0.856 in substantia nigra, using a randomisation approach). Supplementary Data 10). There was no enrichment in Parkin- To investigate the extent to which ASE sites tagged gene-level son’s disease or schizophrenia heritability among eQTLs using the or transcript-specific cis-regulatory effects, we focussed on same method. Thus, although we recognise that eQTLs can be common ASE sites (seen in ≥10 individuals) that were powerful when linked to even more specific cell types for this type unidirectional in nature (same direction of effect across all of analysis20,36, we demonstrate the additional power of ASE individuals). For each valid ASE site, we measured exon analysis to generate disease-relevant information, despite the expression across all three genotypes. Of all testable ASEs, we small number of samples we had at our disposal. found that 43.2% were also likely to be eQTLs. To ask whether the underlying cis-regulatory effects operated in an exon-specific or gene-level manner, we repeated the analysis using gene-level Discussion expression across the genotypes. Of the ASEs that were also likely The human brain is an especially challenging organ in which to eQTLs, 51.8% appeared to operate in an exon-specific manner, conduct eQTL and ASE studies. In addition to the difficulties of implying that they tagged splicing eQTLs. rs7724759, a splice site sample collection, the brain is a highly complex organ, with site- variant present in the CAST gene (Supplementary Fig. 6), and specific pathologies that motivate the use of equivalently specific rs1050078, a variant in SNX19 (Supplementary Fig. 6), are analyses. The brain is also known to express many transcripts not examples of ASEs likely to be driven by exon-level and gene-level seen in other parts of the body, and it is suspected that much of regulation respectively. its transcriptome remains uncharacterised16,37. Although this approach allowed us to identify ASE sites tagging We tackled these issues by collecting RNA-seq data from splicing eQTLs, it was limited to a small subset of common ASE human substantia nigra and putamen, and applied a bank of five sites and represented only 0.8% of all ASEs. To address this issue, transcriptome quantification methods, including annotation- we used the machine learning programme SPIDEX to predict the agnostic approaches as well as approaches to interrogate differ- effect of all ASEs on splicing33. Given a genetic variant, SPIDEX ent stages of RNA processing. We found that there is significant provides the delta percent inclusion ratio (ΔΨ) for the exon in variation among eQTL classes in their neuron- and brain-specific which the variant is located (reported as the maximum ΔΨ across information content, as measured by the cell type-specific tissues). We compared predicted ΔΨ values at ASE sites versus enrichment of eQTL targets. The most neuronally enriched and non-ASE sites and found significantly higher values among ASEs brain-specific results were found in eQTL classes that most NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x a Substantia nigra Putamen Astrocyte Cell type Endothelial Astrocyte Endothelial Cell type Microglia Microglia Neuron Oligodendrocyte Neuron Oligodendrocyte 0 25 50 75 100 0 50 100 150 200 –log10 (p value) –log10 (p value) b Allele-specific expression Substantia nigra Putamen PD SCZ PD SCZ PD SCZ All All Both e-eQTL e-eQTL PUTM g-eQTL g-eQTL SNIG i-eQTL i-eQTL 0 1 2 3 0 1 2 3 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 0.0 0.5 1.0 1.5 2.0 Coefficient p -value (–log10 p value) Coefficient p -value (–log10 p value) Coefficient p -value (–log10 p value) Both PUTM SNIG All e-eQTL g-eQTL i-eQTL Fig. 7 Allele-specific expression sites biologically and disease-relevant information. a Genic locations of ASEs were highly enriched for genes with cell- biased expression in both putamen and substantia nigra. Enrichment of genes with cell-biased expression within ASE locations was performed separately for each tissue and was determined using a Fisher’s Exact test and a significance cutoff of p value < 0.05 (dashed red line at −log10(P) = 1.30). b Enrichment of heritability for Parkinson’s Disease and schizophrenia in ASEs and eQTLs identified in substantia nigra, putamen and across both tissues. Source data are provided as a Source Data file. closely tag the regulation of splicing (e-eQTLs and ex-ex-eQTLs) overlap with enhancer regions, suggesting that many could also rather than gene-level expression. This finding is consistent with represent eRNAs (3.03-fold enrichment over expressed exons). recent studies that suggest that splicing eQTLs can provide sig- Finally, we show that intergenic eQTLs (i-eQTLs) are enriched nificant insights into complex diseases in general38, and brain- for neuronally relevant information, and most importantly that related disorders in particular (e.g., schizophrenia20,36). Thus, in they can provide unique disease insights that would be missed addition to providing a rich eQTL resource, our study suggests using standard analyses, as illustrated by the colocalisation of i- that the utility of existing and future eQTL analyses in human eQTL signals with schizophrenia risk loci. brain may critically depend on the ability of the RNA sequencing Nevertheless, the identification of splicing eQTLs from technology, and of the analytic methods applied, to capture homogenates of macro-dissected human brain, particularly from transcript-specific information. brain regions that are hard to obtain in large numbers, is likely to We also asked whether the incomplete annotation of the remain challenging even after accounting for the on-going human brain transcriptome might be limiting eQTL discovery as development of tools to optimise transcriptome quantification. well as reducing the tissue-specific nature of the regulatory This motivates the use of ASE analysis, a form of within- information discovered. We focused on transcription within individual comparison that compares variation in expression intergenic regions, as transcriptional activity in these regions between two haplotypes of a diploid individual. This within- cannot be explained by the presence of pre-mRNA but instead individual comparison means that ASE analysis is unaffected by could be generated through the expression of long intergenic between-individual confounders, such as the variability in cell non-coding RNAs and enhancer RNAs, which are reported to be type-specific density among individuals. We applied ASE analysis expressed in a highly tissue-specific manner15,39. We show that to 49 putamen and 35 substantia nigra samples, for which both these expressed intergenic regions are reliably detected, and that whole-exome sequencing and genotyping data were available. ~16.1% of these expressed regions are highly likely to represent Consistent with our expectation, we found that ASEs were sig- novel exons of known genes (as demonstrated through the exis- nificantly enriched for variants identified as splicing eQTLs tence of junction spanning reads). They are also enriched for within our own analysis or predicted to affect splicing according 10 NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x ARTICLE to SPIDEX. Furthermore, we found that the ASEs we identified heterozygotes present and variants with mismatching alleles from 1000 Genomes tagged regulatory information that was highly enriched for neu- Project were removed from the analysis. Imputation was performed using MaCH47 and minimac48 using the European panel of the 1000 Genomes Project (March rons and brain-relevant cell types, even after accounting for the 2012: Integrated Phase I haplotype release version 3, based on the 2010 November general enrichment in brain-specific information contained data freeze and 14 March 2012 haplotypes). We used the resulting 5,878,211 SNPs within the RNA-seq data. Finally, and most importantly, we used and 576,942 indels with good postimputation quality (R2 > 0.50) and minor allele two separate approaches to demonstrate the relevance of ASEs to frequency (MAF) of at least 5%. both Parkinson’s disease and schizophrenia, with evidence for A subset of DNA samples were also analysed using whole-exome sequencing49 (N = 57, Supplementary Data 1, EGAS00001002113). Exome sequencing was enriched heritability among ASEs. Given the small numbers of performed on each subject of the cohort according to the manufacturer’s capture samples used in our ASE analyses, this finding is particularly protocol; capture kits included Illumina and Nimblegen. Paired-end sequence reads striking. Thus, we provide evidence to suggest that ASE analysis were aligned with BWA against the reference human genome (UCSC hg19)50. may be a particularly effective and efficient means of obtaining Duplicate read removal, format conversion, and indexing was performed with Picard (http://picard.sourceforge.net/). The Genome Analysis Toolkit was used to regulatory information relevant to splicing, cell type and disease. recalibrate base quality scores, perform local realignments around indels, and to In summary, by using a range of methods to quantify and call and filter the variants51,52. SnpEff was used to annotate gene and effect analyse brain transcriptomic data, we demonstrate the impor- information for the variants31. Subject QC was performed using typical methods to tance of capturing information on the regulation of known and check for call rate, heterozygosity outliers, gender, relatedness, and population outliers. These QC checks were performed using Plink53 based on the Single novel splicing for the understanding of complex brain disorders, Nucleotide Variants (SNVs) that are in the intersection of the exome capture kits and show that ASE analyses performed even on small sample sets used. VCFTools (http://vcftools.sourceforge.net/) was used to convert the SNV vcf can provide additional insights. file to Plink formatted genotypes. Population outlier checks included subjects, populations, and genotypes from the 1000 Genomes Project54, Phase1 Release v2.20101123 for reference. Individual genotypes were removed with genotype Methods quality Phred-scores below 40. Omni-1M, Immunochip and exome sequencing Generation and processing of RNA sequencing data. Human brain samples data were then merged with priority placed on the genotyping array data, only originating from 117 individuals of European descent were obtained from the overwriting calls from arrays when genotyping data were missing. MaCH47 was Medical Research Council (MRC) Sudden Death Brain and Tissue Bank and the used to determine phasing and the resulting files were converted into a variant call Sun Health Research Institute. All samples were authorised for ethically approved format (VCF) file using R scripts. An average of 286,824 SNPs with heterozygous scientific investigation (Research Ethics Committee number 10/H0716/3, National genotypes were obtained per individual. Hospital for Neurology and Neurosurgery and Institute of Neurology Research Ethics Committee) and had fully informed consent for retrieval. These samples constituted a subset of the United Kingdom Brain Expression Consortium eQTL discovery and replication. Prior to performing eQTL analyses, the Prob- (UKBEC) data set. RNA was isolated by using the miRNeasy 96 sample kit abilistic Estimation of Expression Residuals55 (PEER) method was used to identify (Qiagen, UK) and the RNA integrity number was assessed for each sample using unknown factors affecting expression levels and so optimise eQTL discovery. PEER the Agilent 2100 Bioanalyzer (Agilent Technologies UK Ltd, UK) in combination was run using default parameters with brain regions, age, and gender accounted for with the RNA 6000 Nano-LabChip kit (Supplementary Table 1). as known covariates to generate 13 unknown factors (captured through the PEER cDNA libraries were prepared by the UK Brain Expression Consortium in axes), which were applied to RPKM normalised values to produce residuals. As conjunction with AROS Applied Biotechnology A/S (Aarhus, Denmark). Using would be expected, many of the 13 unknown factors were highly correlated with 100 ng of total RNA as input, amplified cDNA was generated with the NuGen recognised covariates such as RIN, library batch, and intronic read rate (Supple- Ovation RNA-seq System V2 (NuGen Technologies, US) according to the mentary Fig. 5). We also identified nine factors that were highly correlated with the manufacturer’s protocol. This protocol deselects for rRNA by using oligo dT and specific brain region analysed. We found that expression of SLC6A3, which encodes random hexamer primers for reverse transcription. After fragmentation of cDNA the dopamine transporter (expressed with very high specificity in dopaminergic (1 μg) using a Covaris S220 Ultrasonicator, we used the Illumina TruSeq DNA neurons of the ventral midbrain including the substantia nigra, see below), was library preparation kit (Illumina, US) followed by 10 cycles of PCR amplification of highly correlated with the 10th PEER axis (p value: 1.9 × 10−26). library molecules containing adaptor molecules on both ends to generate the final Variants within ±1 Mb of each normalised expression feature (i.e., gene, exon, libraries. Sequencing of these DNA libraries was performed with Illumina’s TruSeq exon–exon junction, and unannotated intergenic region) were tested for cis-eQTL V3 chemistry/HiSeq2000 to generate 100 base pair paired-end reads. Illumina’s discovery using the R package MatrixEQTL56. Gender, age, and the first three CASAVA Software was then used to generate fastq-files. Paired-end data were genetic principal component vectors were included as covariates within the linear mapped to the human genome (build GRCh37) using tophat240 (v2.0.9) with model of marker genotype (imputed expected counts of minor allele) against PEER default settings and a transcriptome-guided approach using the Ensembl reference corrected expression values. The Benjamini–Hochberg method was applied to (v72) based on GENCODE version 18. Reads mapping to rRNA regions were calculate the FDR adjusting for the tests performed for each transcriptomic feature. removed from the analysis. Finally, we used a stepwise conditional analysis for each quantification type to Transcriptome quantification was performed using the transcriptome definition detect independent variant effects targeting the same expression feature. from Ensembl reference (v72) using HTSeq-Counts41 (v0.5.4p4) for exonic regions, Three major types of data were used to validate eQTL results: (i) eQTL data BEDtools42 for intronic regions, DEXSeq43 (v1.10.6) for individual exons, Altrans44 reported by Ramasamy and colleagues17, which used an overlapping set of donor (v1.1.02) for exon–exon junctions and the derfinder R package45 for unannotated samples but a hybridisation-based microarray method for exon-level transcriptome transcribed regions. Expressed regions were identified by derfinder using default quantification, (ii) eQTL data generated by the GTEx19, PsychENCODE and settings and filtering for expressed regions of >100 bp in length, annotated as CommonMind Consortia20,21, which were based on independent sample sets, but intergenic based on Ensembl v72 and UCSC through the R library TxDb.Hsapiens. assayed multiple brain regions and which used RNA sequencing to quantify the UCSC.hg19.knownGene version v3.1.2, and which uniquely mapped (defined as transcriptome, and (iii) eQTL data relating to 373 lymphoblastoid cell lines >98% alignment precision). This approach was applied to each tissue separately generated by the Geuvadis consortium10, which used RNA sequencing for (putamen and substantia nigra). All forms of transcriptome quantification were transcriptome quantification. Gene-level eQTLs (ge-eQTLs) identified within our normalised using CQN46, with GC content calculated separately for each type of study were declared as replicated when the significant SNP-gene pair was declared quantification and used as an input for the CQN R-bioconductor package. The a hit with an FDR < 5% in both our data set and the comparison data set. resulting normalised expression data were then transformed into Reads Per Kilobase of transcript per Million (RPKM) values and log2 converted. Exonic, Characterisation of unannotated transcribed intergenic regions. We leveraged intronic, and transcribed intergenic regions with an RPKM of >0.1 in at least 80% information from split reads (reads aligning to a genomic location with a gap or of samples were selected for downstream analyses. multiple gaps) and combined this with information on co-expression and physical proximity to known genes to characterise novel transcribed regions. In the first DNA extraction and genotyping. DNA was extracted from sub-dissected samples instance, the split read information for each sample was collected from the (100–200 mg) of human post mortem brain tissue using Qiagen’s DNeasy Blood & Tophat240 junction output file, and overlaps between transcribed intergenic regions Tissue Kit (Qiagen, UK). All samples were genotyped on the Illumina Infinium targeted by eQTLs and split reads were assessed using the GenomicRanges57 R Omni1-Quad BeadChip and on the Immunochip, a custom genotyping array package. If a transcribed region overlapped with a split read, then the same split designed for the fine mapping of autoimmune disorders. Standard quality control read was screened for overlap with known genes. Only transcribed intergenic on the merged genotyping data was performed. Individuals of suspected non- regions with split reads present in at least four separate samples were classed as European descent and samples with percentage of non-missing genotypes of < 95% having high evidence for being part of a known gene. In the absence of relevant were removed from the analysis. Reported gender status and non relatedness of split read data, the nearest gene (in genomic distance) from each transcribed samples were confirmed. Monomorphic SNPs, variants with missing position intergenic region was selected using the GenomicRanges57 R package and the r² information, variants with a p value < 0.0001 for deviation from Hardy–Weinberg was calculated between the expression of the transcribed intergenic region and each equilibrium, variants with a genotype call rate < 95%, variants with less than two exon of the nearest gene. Transcribed regions that had a maximum r² > 0.2 and NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x were <5Kb from the nearest gene were classed as having moderate evidence for ordered alphabetically, were used to calculate the adjusted ratio (C1 + 0.5)/(C2 + being part of a known gene. Transcribed intergenic regions were classed as having 0.5) to measure the ASE effect size. The use of +0.5 as an adjustment was moti- low evidence for being part of a known gene when the r2 was <0.2 or they were vated by the predicted reduction in the Taylor-Maclaurin bias estimator. Statistical >5Kb from the nearest gene. tests for ASE were carried out via exact tests for departure from binomial expec- We also explored the possibility that unannotated transcribed regions could tation under the null hypothesis of (C1 = C2), using the binom.test() function in R. represent enhancer RNAs (eRNAs). We used the GeneHancer database v4.423, Tests were only conducted at heterozygous SNP sites where the total number of which draws information from a wide range of sources (including ENCODE, the counts (C1 + C2) exceeded five reads (defining valid sites). For each sample, we FANTOM5 atlas, the VISTA Enhance Browser, dbSUPER and EPDnew, and converted the p values into FDR using the Benjamini–Hochberg procedure. An UCNEbase) to define enhancer locations and then quantified the percentage FDR threshold of 5% was used to assess ASE significance in the valid sites. ASE overlap between eQTL target regions and enhancer regions. We adjusted for signal validation was performed using ASE data generated from previously pub- differences in the genomic size of eQTL target regions by calculating the percentage lished lymphoblastoid cell line data10. After filtering for SNPs analysed in both data overlap per Mb of transcribed target sequence. sets (n = 54,214) we declared as validated those ASEs discovered in our brain data set that also had an FDR-corrected p value of <5% within the lymphoblastoid data set. Validation of unannotated transcribed intergenic regions. We validated the transcription of intergenic regions targeted by i-eQTLs in silico using the data available on tissue-specific transcription generated by the GTEx consortium (www. Characterisation of ASE signals. Variant Effect Predictor (VEP)62, SNPeff31, and gtexportal.org) and mapped through recount222. The genome coordinates of all SPIDEX33 were used to annotate all hetSNPs and to obtain a list of genes that had intergenic regions of interest were quantified from the transcription expression at least one significant hetSNP. For a subset of ASEs, namely those annotated as profiles generated with derfinder and available in recount222. Counts were trans- protein-truncating variants within VEP (defined as stop gain, donor splice site and formed to RPKMs and validation of transcription was considered when the region donor acceptor mutations) or common across individuals (defined as a ASEs had an RPKM > 0.1 in at least 80% of samples in the tissue of interest. identified in ≥10 individuals), we also investigated the impact of genotype on gene We selected eight transcribed intergenic regions targeted by i-eQTLs for expression. The average expression of the exon or gene (following PEER correc- validation by Sanger sequencing. All regions were detected through analysis of tion) was calculated separately for individuals homozygous or heterozygous for the putamen samples and we used a subset of RNA samples from the set of 111 used variant of interest or if possible across all three genotypes. We then performed for RNA sequencing (Supplementary Table 2). In each case, reverse transcription t tests or ANOVA as appropriate to identify departures from the null hypothesis of was performed with 500 μg of total RNA using High Capacity cDNA RT Kit no association between expression and genotype. (Applied Biosystems) and random primers as per manufacturer’s instructions. PCR was performed using specific primers (Supplementary Table 2), all of which were designed to span predicted exon–exon junctions, and FastStart PCR Master Investigation of the disease relevance of ASEs using randomisation. We (Roche). Amplification of the predicted band was confirmed by agarose gel investigated the enrichment of GWAS risk loci for Parkinson’s disease (with the electrophoresis and following confirmation, enzymatic clean-up of PCR products exclusion of data generated by 23andMe)1 and schizophrenia7 among ASEs versus was performed using Exonuclease I (Thermo Scientific) and FastAP non-ASE sites by using a randomisation approach similar to that described by Thermosensitive Alkaline Phosphatase (Thermo Scientific). Sequencing was Nicolae and colleagues34 to generate an empirical p value for the significance of performed using the BigDye terminator kit (Applied Biosystems). All sequences overlaps between risk SNPs and ASE sites, accounting for differences in read depth. were viewed using CodonCode Aligner (V. 6.0.2). Let C be the set of SNPs that are in common between our set of hetSNPs (used for our ASE analyses) and the set of variants available from the GWAS summary statistics of the disease association study. Let Ase be the subset of SNPs in C Identification of redundant eQTLs by beta-hetereogeneity testing. A mixed declared to be ASE hits in our study, and let NonAse be the subset of SNPs in C not model approach was used to identify heterogeneity in eQTL signal strength (beta declared to be ASE hits. For the randomisation procedure, we repeated the coefficient or slope) within a gene. To test for beta-heterogeneity between ge-eQTL following two steps 105 times: and each other eQTL class (namely gi-eQTLs, e-eQTLs, ex-ex-eQTLs and relevant i-eQTLs) targeting the same gene, two models were fitted: (1) a non-heterogeneous i. Randomly select SNPs with replacement from NonAse to create a SNP set (single slope) model with allele dosage as the main effect and two random effects, with the same size and read depth distribution as that of Ase. This was done one indexing individuals and one indexing genes or exons or exon junctions and by placing the NonAse into bins based on their average read depth across (2) a heterogeneous (multiple slope) model containing the same terms, but with the samples and randomly selecting SNPs from the same to match the average addition of a fixed-effect allele dosage × exon-or-exon junction index interaction read distribution of the Ase set. term. The lme4 R package was used to fit both models. A p value for the likelihood ii. Calculate the enrichment of GWAS signals within this SNP set by recording ratio test comparing the two models was generated with the R function anova() and the GWAS p values for each SNP within the set and computing the mean of an FDR of <5% was applied to assess significance. the −log10 transformed values. Let ri be the mean value for the ith iteration. The z score for Ase (zAse) was obtained using the mean and standard deviation Identification of disease-relevant eQTLs. We used the STOPGAP database25 of the ri values, and was tested using the pnorm() function in R to obtain a p value (accessed on 20th of March 2018) to access and sub-classify loci identified through for the enrichment of GWAS risk SNPs among ASEs versus non-ASEs. genome-wide association studies (GWAS). eQTL-GWAS overlap was checked Using a similar approach, we also tested ASEs for evidence of enrichment of using all eQTLs passing an FDR < 5%. The percentage overlap and enrichment was eQTLs and SNPs with effects on exon inclusion as predicted by SPIDEX. calculated based on the total number of eQTLs identified. In addition, summary statistics were obtained from Parkinson’s disease (with the exclusion of data gen- Assessment of brain cell type specificity for eQTLs and ASEs. We used the erated by 23andMe) and schizophrenia GWAS1,7. We applied coloc26 to colocalise WGCNA R package63 to construct separate gene co-expression networks for Parkinson’s disease and schizophrenia loci with our eQTL signals. For each locus substantia nigra and putamen in order to generate modules enriched for cell type- with a GWAS p value of <1 × 10−5 we examined all SNPs available in both data sets specific gene expression profiles. As input for each network construction, we used within 1 Mb of the GWAS SNP of interest, and ran coloc with default parameters genes detected in >70% of samples and with an RPKM of >0.1. Following nor- and priors. We called the signals colocalised when Coloc H3 + H4 probabilities malisation using CQN46 and accounting for GC content, we used 13 PEER axes were greater than 0.8 and the H3/H4 probability ratio was greater than 2. (see above) in addition to gender, age, and genetic axes to perform covariate correction. The number of PEER axes was optimised to maximise the clustering of ASE signal discovery and replication. The ASE discovery pipeline was motivated high quality cell markers64. We generated signed networks using a soft- by the method of Rozowsky and colleagues58 and involved the creation of parent thresholding power of 7 for substantia nigra and 8 for putamen networks, which haploids using phased variant data to reduce the impact of mapping biases on ASE approximated a scale-free topology. We obtained first-pass definitions of the identification. Each individual’s haploid genome was constructed using genotype modules via the dynamic Tree Cutting algorithm in WGCNA, and we then refined data in the relevant VCF file, as parental information was unavailable. Deviations these modules by applying a k-means algorithm as described in our previous from the reference genome present in the VCF file were used to update the publication65. Modules were annotated in terms of cell-specific enrichments using reference and create an artificial genome (two haploid genomes originating from the userListEnrichment R function implemented in the WGCNA R Package, each parent) using the Personal Genome Constructor tool vcf2diploid58 (version combined with additional enrichment analysis on other marker gene sets64,66–68. 0.2.4). The haploid genomes were arbitrarily referred to as parent1 and parent2. Consistent with other studies, we robustly identified modules enriched for cell Trimmed fastq reads were aligned to individual haploid parent genomes using type-specific gene expression profiles21,66–68. Genes assigned to modules sig- Tophat40 and Bowtie259 (version 2.0.6), following a transcriptome-guided nificantly enriched for brain-related cell type markers and with a module mem- approach, setting parameters to acquire the best alignment and allowing up to two bership of >0.3 were allocated a cell type label of neuron, microglia, astrocyte, mismatches. The alignments acquired for both haploid genomes, for each indivi- oligodendrocyte, and endothelial. dual, were merged using the Suspenders tool (version 0.2.3), selecting the single Next, for each eQTL targeting a known genic region, if the target gene was best-quality alignment. IGVTools60,61, version 2.3.18, was used to count the reads allocated to a cell type then the related eQTL received the same cell type label. In with one or the other hetSNP allele. Counts of both alleles, C1 and C2, with alleles the case of eQTLs targeting unannotated transcribed intergenic regions with high 12 NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x ARTICLE or moderate evidence for association to a known gene, the eQTL received a cell 11. Kim-Hellmuth, S. et al. Genetic regulatory effects modified by immune type label based on that known gene’s network position. For eQTLs targeting activation contribute to autoimmune disease associations. Nat. Commun. 8, transcribed intergenic regions with low evidence for association with a known gene 266 (2017). or which could not be classified, we assigned the target expression feature to a 12. Kochi, Y. Genetics of autoimmune diseases: perspectives from genome-wide module (and by inference a cell type) based on its highest module membership association studies: Table 1. Int. Immunol. 28, 155–161 (2016). (defined as the correlation of expression with the first principal component 13. Raj, T. et al. Polarization of the effects of autoimmune and neurodegenerative (“eigengene”) of each module), provided the module membership was at least 0.3. risk alleles in leukocytes. Science 344, 519–523 (2014). Finally, for each eQTL class and each cell type we applied a Fisher’s Exact test to 14. Hawrylycz, M. J. et al. An anatomically comprehensive atlas of the adult test for enrichment of that cell type label among the genes associated to the eQTL human brain transcriptome. Nature 489, 391–399 (2012). class, relative to the genes not associated with that eQTL class but with module 15. Briggs, J. A., Wolvetang, E. J., Mattick, J. S., Rinn, J. L. & Barry, G. membership >0.3 to a relevant cell type module. We followed a similar approach Mechanisms of long non-coding rnas in mammalian nervous system for ASEs and assessed ASE-containing genes for significant cell type enrichment. development, plasticity, disease, and evolution. Neuron 88, 861–877 (2015). 16. Jaffe, A. E. et al. Developmental regulation of human cortex transcription and Partitioned heritability analysis of ASE sites. GWAS summary statistics were its clinical relevance at single base resolution. Nat. Neurosci. 18, 154–161 obtained for schizophrenia7 and Parkinson’s disease (with the exclusion of data (2015). generated by 23andMe)1. Stratified LD score regression67, a method for parti- 17. Ramasamy, A. et al. Genetic variability in the regulation of gene expression in tioning SNP heritability across functional genomic annotations, was used to test for ten regions of the human brain. Nat. Neurosci. 17, 1418–1428 (2014). enrichment of schizophrenia and Parkinson’s disease heritability within ASE 18. Kang, H. J. et al. Spatio-temporal transcriptome of the human brain. Nature annotations. Region-specific annotations were defined for putamen and substantia 478, 483–489 (2011). nigra (ASEs with FDR < 0.05 within a brain region were assigned to that region). In 19. Ardlie, K. G. et al. The Genotype-Tissue Expression (GTEx) pilot analysis: addition, an annotation (all) containing all ASEs irrespective of brain region with multitissue gene regulation in humans. Science 348, 648–660 (2015). FDR < 0.05 was defined. Annotations were added individually to the “baseline” 20. Fromer, M. et al. Gene expression elucidates functional impact of polygenic model of 53 annotations provided by Finucane et al., which comprises 24 different risk for schizophrenia. Nat. Neurosci. 19, 1442–1453 (2016). genome-wide annotations reflecting genetic architecture, such as conserved regions 21. Wang, D. et al. Comprehensive functional genomic resource and integrative and histone marks. HapMap Project Phase 3 SNPs and 1000 Genome Project model for the human brain. Science 362, eaat8464 (2018). European population SNPs were used for the regression and LD reference panel, 22. Collado-Torres, L. et al. Reproducible RNA-seq analysis using recount2. Nat. respectively. Only SNPs with MAF > 5% were used for heritability partitioning and Biotechnol. 35, 319–321 (2017). the HLA region was excluded. We report the regression coefficient p value, which 23. Fishilevich, S. et al. GeneHancer: genome-wide integration of enhancers and tests whether the ASE annotations contribute significantly to SNP heritability after target genes in GeneCards. Database (Oxford). 2017 bax028 (2017). controlling for the effects of the “baseline” model. 24. Jaffe, A. E. et al. Developmental and genetic regulation of the human cortex transcriptome illuminate schizophrenia pathogenesis. Nat. Neurosci. 21, Reporting summary. Further information on research design is available in 1117–1125 (2018). the Nature Research Reporting Summary linked to this article. 25. Shen, J., Song, K., Slater, A. J., Ferrero, E. & Nelson, M. R. STOPGAP: a database for systematic target opportunity assessment by genetic association predictions. Bioinformatics 33, 2784–2786 (2017). Data availability 26. Giambartolomei, C. et al. Bayesian test for colocalisation between pairs of We have made data available in a web server and via full and summary-statistic data genetic association studies using summary statistics. PLoS Genet. 10, e1004383 downloads. Using our web resource, http://braineacv2.inf.um.es/, users can access and (2014). visualise all forms of transcriptome quantification, eQTLs as well as gene co-expression 27. Baran, Y. et al. The landscape of genomic imprinting across diverse adult networks (Supplementary Fig. 1). The RNA-seq, whole-exome sequencing and human tissues. Genome Res. 25, 927–936 (2015). genotyping data can be accessed through the European Genome-phenome Archive 28. Docherty, L. E. et al. Genome-wide DNA methylation analysis of patients with numbers EGAS00001002113 and EGAS00001003065. The source data underlying imprinting disorders identifies differentially methylated regions associated Figs. 1b, 2b–d, 3a–d, 4b, 5a, b, 6b, c, 7b, and 8a and Supplementary Figs. 2a–c and 3–6 are with novel candidate imprinted genes. J. Med. Genet. 51, 229–238 (2014). provided as a Source Data file. 29. Monk, D. et al. Human imprinted retrogenes exhibit non-canonical imprint chromatin signatures and reside in non-imprinted host genes. Nucleic Acids Received: 5 April 2019; Accepted: 20 December 2019; Res. 39, 4577–4586 (2011). 30. Ning, Z. et al. Regulation of SPRY3 by X chromosome and PAR2-linked promoters in an autism susceptibility region. Hum. Mol. Genet. 24, 5126–5141 (2015). 31. Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly. (Austin). 6, 80–92 (2012). References 32. Castel, S. E., Levy-Moonshine, A., Mohammadi, P., Banks, E. & Lappalainen, 1. Chang, D. et al. A meta-analysis of genome-wide association studies identifies T. Tools and best practices for data processing in allelic expression analysis. 17 new Parkinson’s disease risk loci. Nat. Genet. 49, 1511–1516 (2017). Genome Biol. 16, 195 (2015). 2. Nalls, M. A. et al. Imputation of sequence variants for identification of genetic risks for Parkinson’s disease: a meta-analysis of genome-wide association 33. Xiong, H. Y. et al. RNA splicing. The human splicing code reveals new insights into the genetic determinants of disease. Science 347, 1254806 (2015). studies. Lancet 377, 641–649 (2011). 34. Nicolae, D. L. et al. Trait-associated SNPs are more likely to be eQTLs: 3. McKenzie, M. et al. Overlap of expression Quantitative Trait Loci (eQTL) in annotation to enhance discovery from GWAS. PLoS Genet. 6, e1000888 human brain and blood. BMC Med. Genomics 7, 31 (2014). 4. Harold, D. et al. Genome-wide association study identifies variants at CLU (2010). 35. Finucane, H. K. et al. Partitioning heritability by functional annotation using and PICALM associated with Alzheimer’s disease. Nat. Genet. 41, 1088–1093 genome-wide association summary statistics. Nat. Genet. 47, 1228–1235 (2009). (2015). 5. Hollingworth, P. et al. Common variants at ABCA7, MS4A6A/MS4A4E, EPHA1, CD33 and CD2AP are associated with Alzheimer’s disease. Nat. 36. Takata, A., Matsumoto, N. & Kato, T. Genome-wide identification of splicing QTLs in the human brain and their enrichment among schizophrenia- Genet. 43, 429–435 (2011). associated loci. Nat. Commun. 8, 14519 (2017). 6. Lambert, J.-C. et al. Genome-wide association study identifies variants at CLU 37. Zhang, D. et al. Incomplete annotation of OMIM genes is likely to be limiting and CR1 associated with Alzheimer’s disease. Nat. Genet. 41, 1094–1099 the diagnostic yield of genetic testing, particularly for neurogenetic disorders. (2009). Preprint at https://doi.org/10.1101/499103v2 (2018). 7. Pardiñas, A. F. et al. Common schizophrenia alleles are enriched in mutation- 38. Li, Y. I. et al. RNA splicing is a primary link between genetic variation and intolerant genes and in regions under strong background selection. Nat. Genet. disease. Science 352, 600–604 (2016). 50, 381–389 (2018). 39. Yao, P. et al. Coexpression networks identify brain region-specific enhancer 8. Ripke, S. et al. Biological insights from 108 schizophrenia-associated genetic loci. Nature 511, 421–427 (2014). RNAs in the human brain. Nat. Neurosci. 18, 1168–1174 (2015). 40. Kim, D. et al. TopHat2: accurate alignment of transcriptomes in the presence 9. Westra, H.-J. & Franke, L. From genome to function by studying eQTLs. of insertions, deletions and gene fusions. Genome Biol. 14, R36 (2013). Biochim. Biophys. Acta—Mol. Basis Dis. 1842, 1896–1902 (2014). 41. Anders, S., Pyl, P. T. & Huber, W. Genome analysis HTSeq - a Python 10. Lappalainen, T. et al. Transcriptome and genome sequencing uncovers functional variation in humans. Nature 501, 506–511 (2013). framework to work with high-throughput sequencing data. 31, 166–169 (2014). NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications 13 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x 42. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for Acknowledgements comparing genomic features. Bioinformatics 26, 841–842 (2010). Mina Ryten, David Zhang, and Karishma D’Sa were supported by the UK Medical Research 43. Anders, S., Reyes, A. & Huber, W. Detecting differential usage of exons from Council (MRC) through the award of Tenure-track Clinician Scientist Fellowship to Mina RNA-seq data. Genome Res. 22, 2008–2017 (2012). Ryten (MR/N008324/1). Sebastian Guelfi was supported by Alzheimer’s Research UK 44. Ongen, H. & Dermitzakis, E. T. Alternative splicing QTLs in European and through the award of a PhD Fellowship (ARUK-PhD2014-16). Regina Reynolds was sup- African populations. Am. J. Hum. Genet. 97, 567–575 (2015). ported through the award of a Leonard Wolfson Doctoral Training Fellowship in Neuro- 45. Collado-Torres, L. et al. Flexible expressed region analysis for RNA-seq with degeneration. All RNA sequencing data performed as part of this study were generated by derfinder. Nucleic Acids Res. 45, e9–e9 (2017). the commercial company AROS Applied Biotechnology A/S (Denmark). 46. Hansen, K. D., Irizarry, R. A. & Wu, Z. Removing technical variability in RNA-seq data using conditional quantile normalization. Biostatistics 13, 204–216 (2012). Author contributions Sebastian Guelfi processed the RNA sequencing data and generated the eQTL data. 47. Li, Y., Willer, C. J., Ding, J., Scheet, P. & Abecasis, G. R. MaCH: using Karishma D’Sa generated the ASE data and with Juan Botia developed and implemented sequence and genotype data to estimate haplotypes and unobserved the randomisation analysis. Juan Botia also performed all network-based analyses. Jana genotypes. Genet. Epidemiol. 34, 816–834 (2010). Vandracova assisted with the annotation of eQTLs and ASEs. David Zhang performed 48. Howie, B., Fuchsberger, C., Stephens, M., Marchini, J. & Abecasis, G. R. Fast coloc analyses. Regina Reynolds validated novel transcribed regions and with the assis- and accurate genotype imputation in genome-wide association studies tance of Sarah Gagliano conducted analyses using LD score regression. Adaikalavan through pre-phasing. Nat. Genet. 44, 955–959 (2012). Ramasamy performed the analysis of the genotyping data and contributed to the 49. Jansen, I. E. et al. Discovery and functional prioritization of Parkinson’s development of the ASE and eQTL pipelines. Daniah Trabzuni performed the RNA disease candidate genes from large-scale whole exome sequencing. Genome extractions from post mortem human brain tissue. Leonardo Collado-Torres assisted Biol. 18, 22 (2017). with the analysis of data available through recount2. Andy Thomason and Pedro Quijada 50. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows- Leyton designed and generated the webserver to enable access to eQTL and ASE data Wheeler transform. Bioinformatics 25, 1754–1760 (2009). generated within this study. Mike Nalls provided assistance with the use of the Par- 51. McKenna, A. et al. The Genome Analysis Toolkit: a MapReduce framework kinson’s Disease GWAS data, which was supplied by the IPDGC. Members of the for analyzing next-generation DNA sequencing data. Genome Res. 20, IPDGC assisted with the preparation of the manuscript. Kerrin Small provided input 1297–1303 (2010). regarding the development of the randomisation analysis. Colin Smith assisted with the 52. DePristo, M. A. et al. A framework for variation discovery and genotyping collection of human putamen and substantia nigra samples together with the UK Brain using next-generation DNA sequencing data. Nat. Genet. 43, 491–498 Expression Consortium. Michael Weale and Mina Ryten supervised the generation of the (2011). ASE and eQTL pipelines. Sebastian Guelfi, Karishma D’Sa, Juan Botia, Jana Vandracova, 53. Purcell, S. et al. PLINK: a tool set for whole-genome association and Michael Weale, John Hardy, and Mina Ryten prepared the manuscript. Michael Weale, population-based linkage analyses. Am. J. Hum. Genet. 81, 559–575 (2007). John Hardy, and Mina Ryten conceived and designed the project. 54. An integrated map of genetic variation from 1,092 human genomes. Nature 491, 56–65 (2012). 55. Stegle, O., Parts, L., Durbin, R. & Winn, J. A bayesian framework to account Competing interests for complex non-genetic factors in gene expression levels greatly increases Author M.E.W. is an employee of Genomics plc, a genomics based healthcare company. power in eQTL studies. PLoS Comput. Biol. 6, 1–11 (2010). His involvement in the conduct of this research was solely in his former capacity as a 56. Shabalin, A. A. Matrix eQTL: ultra fast eQTL analysis via large matrix Reader in Statistical Genetics at King’s College London. operations. Bioinformatics 28, 1353–1358 (2012). 57. Lawrence, M. et al. Software for computing and annotating genomic ranges. PLoS Comput. Biol. 9, e1003118 (2013). Additional information 58. Rozowsky, J. et al. AlleleSeq: analysis of allele-specific expression and binding Supplementary information is available for this paper at https://doi.org/10.1038/s41467- in a network framework. Mol. Syst. Biol. 7, 522 (2011). 020-14483-x. 59. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359 (2012). Correspondence and requests for materials should be addressed to M.R. 60. Robinson, J. T., Thorvaldsdóttir, H., Wenger, A. M., Zehir, A. & Mesirov, J. P. Variant review with the integrative genomics viewer. Cancer Res. 77, e31–e34 Peer review information Nature Communications thanks Ezgi Hacisuleyman and the (2017). other, anonymous, reviewer(s) for their contribution to the peer review of this work. 61. Thorvaldsdottir, H., Robinson, J. T. & Mesirov, J. P. Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Reprints and permission information is available at http://www.nature.com/reprints Brief. Bioinform. 14, 178–192 (2013). 62. McLaren, W. et al. The Ensembl variant effect predictor. Genome Biol. 17, 122 Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in (2016). published maps and institutional affiliations. 63. Langfelder, P. & Horvath, S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9, 559 (2008). 64. Cahoy, J. D. et al. A transcriptome database for astrocytes, neurons, and Open Access This article is licensed under a Creative Commons oligodendrocytes: a new resource for understanding brain development and Attribution 4.0 International License, which permits use, sharing, function. J. Neurosci. 28, 264–278 (2008). adaptation, distribution and reproduction in any medium or format, as long as you give 65. Botía, J. A. et al. An additional k-means clustering step improves the biological appropriate credit to the original author(s) and the source, provide a link to the Creative features of WGCNA gene co-expression networks. BMC Syst. Biol. 11, 47 Commons license, and indicate if changes were made. The images or other third party (2017). material in this article are included in the article’s Creative Commons license, unless 66. Oldham, M. C. et al. Functional organization of the transcriptome in human indicated otherwise in a credit line to the material. If material is not included in the brain. Nat. Neurosci. 11, 1271–1282 (2008). article’s Creative Commons license and your intended use is not permitted by statutory 67. Winden, K. D. et al. The organization of the transcriptional network in regulation or exceeds the permitted use, you will need to obtain permission directly from specific neuronal classes. Mol. Syst. Biol. 5, 291 (2009). the copyright holder. To view a copy of this license, visit http://creativecommons.org/ 68. Lein, E. S. et al. Genome-wide atlas of gene expression in the adult mouse licenses/by/4.0/. brain. Nature 445, 168–176 (2007). 69. Sibley, C. R. et al. Recursive splicing in long vertebrate genes. Nature 521, 371–375 (2015). © The Author(s) 2020 14 NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x ARTICLE International Parkinson’s Disease Genomics Consortium (IPDGC) Alastair J. Noyce1,14, Aude Nicolas8, Mark R. Cookson8, Sara Bandres-Ciga8, J. Raphael Gibbs8, Dena G. Hernandez8, Andrew B. Singleton8, Xylena Reed8, Hampton Leonard8, Cornelis Blauwendraat8,14, Faraz Faghri3,15,16, Jose Bras10, Rita Guerreiro17, Arianna Tucci17, Demis A. Kia17, Henry Houlden17, Helene Plun-Favreau17, Kin Y Mok17, Nicholas W. Wood17, Ruth Lovering17, Lea R’Bibo17, Mie Rizig17, Viorica Chelban17, Manuela Tan12, Huw R. Morris18, Ben Middlehurst19, John Quinn19, Kimberley Billingsley19, Peter Holmans20, Kerri J. Kinghorn21, Patrick Lewis22, Valentina Escott-Price23, Nigel Williams23, Thomas Foltynie24, Alexis Brice25, Fabrice Danjou25, Suzanne Lesage25, Jean-Christophe Corvol25, Maria Martinez25,26, Anamika Giri27,28, Claudia Schulte27,28, Kathrin Brockmann27,28, Javier Simón-Sánchez27,28, Peter Heutink27,28, Thomas Gasser27,28, Patrizia Rizzu28, Manu Sharma29, Joshua M. Shulman30,31, Laurie Robak30, Steven Lubbe32, Niccolo E. Mencacci33, Steven Finkbeiner34,35, Codrin Lungu36, Sonja W. Scholz37, Ziv Gan-Or38, Guy A. Rouleau39, Lynne Krohan40, Jacobus J. van Hilten41, Johan Marinus41, Astrid D. Adarmes-Gómez42, Inmaculada Bernal-Bernal42, Marta Bonilla-Toribio42, Dolores Buiza-Rueda42, Fátima Carrillo42, Mario Carrión-Claro42, Pablo Mir42, Pilar Gómez-Garre42, Silvia Jesús42, Miguel A. Labrador-Espinosa42, Daniel Macias42, Laura Vargas-González42, Carlota Méndez-del-Barrio42, Teresa Periñán-Tocino42, Cristina Tejera-Parrado42, Monica Diez-Fairen43, Miquel Aguilar43, Ignacio Alvarez43, María Teresa Boungiorno43, Maria Carcel43, Pau Pastor43, Juan Pablo Tartari43, Victoria Alvarez44, Manuel Menéndez González44, Marta Blazquez44, Ciara Garcia44, Esther Suarez-Sanmartin44, Francisco Javier Barrero45, Elisabet Mondragon Rezola46, Jesús Alberto Bergareche Yarza47, Ana Gorostidi Pagola47, Adolfo López de Munain Arregui47, Javier Ruiz-Martínez47, Debora Cerdan48, Jacinto Duarte49, Jordi Clarimón50,51, Oriol Dols-Icardo50,51, Jon Infante51,52,53, Juan Marín51,54, Jaime Kulisevsky51,54, Javier Pagonabarraga51,54, Isabel Gonzalez-Aramburu52, Antonio Sanchez Rodriguez52, María Sierra52, Raquel Duran55, Clara Ruz55, Francisco Vives55, Francisco Escamilla-Sevilla56, Adolfo Mínguez56, Ana Cámara57, Yaroslau Compta57, Mario Ezquerra57, Maria Jose Marti57, Manel Fernández57, Esteban Muñoz57, Rubén Fernández-Santiago57, Eduard Tolosa57, Francesc Valldeoriola57, Pedro García-Ruiz58, Maria Jose Gomez Heredia59, Francisco Perez Errazquin59, Janet Hoenicka60, Adriano Jimenez-Escrig60, Juan Carlos Martínez-Castrillo60, Jose Luis Lopez-Sendon60, Irene Martínez Torres61, Cesar Tabernero61, Lydia Vela61, Alexander Zimprich62, Lasse Pihlstrom63, Sulev Koks64,65,66, Pille Taba67, Kari Majamaa68,69, Ari Siitonen68,69, Njideka U. Okubadejo70 & Oluwadamilola O. Ojo70 14 Preventive Neurology Unit, Wolfson Institute of Preventive Medicine, QMUL, London, UK. 15National Institute of Neurological Disorders and Stroke, Bethesda, MD, USA. 16Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, IL, USA. 17Department of Molecular Neuroscience, UCL, London, UK. 18Department of Clinical Neuroscience, University College London, London, UK. 19Institute of Translational Medicine, University of Liverpool, Liverpool, UK. 20Biostatistics & Bioinformatics Unit, Institute of Psychological Medicine and Clinical Neuroscience, MRC Centre for Neuropsychiatric Genetics & Genomics, Cardiff, UK. 21Institute of Healthy Ageing, University College London, London, UK. 22University of Reading, Reading, UK. 23MRC Centre for Neuropsychiatric Genetics and Genomics, Cardiff University School of Medicine, Cardiff, UK. 24UCL Institute of Neurology, London, UK. 25Institut du Cerveau et de la Moelle épinière, ICM, Inserm U 1127, CNRS, UMR 7225, Sorbonne Universités, UPMC University Paris 06, UMR S 1127, AP-HP, Pitié-Salpêtrière Hospital, Paris, France. 26Paul Sabatier University, Toulouse, France. 27Department for Neurodegenerative Diseases, Hertie Institute for Clinical Brain Research, University of Tübingen, Tübingen, Germany. 28DZNE, German Center for Neurodegenerative Diseases, Tübingen, Germany. 29Centre for Genetic Epidemiology, Institute for Clinical Epidemiology and Applied Biometry, University of Tubingen, Tübingen, Germany. 30Departments of Neurology, Neuroscience, and Molecular & Human Genetics, Baylor College of Medicine, Houston, TX, USA. 31Jan and Dan Duncan Neurological Research Institute, Texas Children’s Hospital, Houston, TX, USA. 32Ken and Ruth Davee Department of Neurology, Northwestern University Feinberg School of Medicine, Chicago, IL, USA. 33Northwestern University Feinberg School of Medicine, Chicago, IL, USA. 34Departments of Neurology and Physiology, University of California, San Francisco, CA, USA. 35Gladstone Institute of Neurological Disease; Taube/Koret Center for Neurodegenerative Disease Research, San Francisco, CA, USA. 36National Institutes of Health Division of Clinical Research, NINDS, National Institutes of Health, Bethesda, MD, USA. 37 Neurodegenerative Diseases Research Unit, National Institute of Neurological Disorders and Stroke, Bethesda, MD, USA. 38Montreal Neurological Institute and Hospital, Department of Neurology & Neurosurgery, Department of Human Genetics, McGill University, Montréal, QC H3A 0G4, Canada. 39Department of Human Genetics, McGill University, Montréal, QC H3A 0G4, Canada. 40Department of Neurology, Leiden University Medical Center, Leiden, Netherlands. 41Instituto de Biomedicina de Sevilla (IBiS), Hospital Universitario Virgen del Rocío/CSIC/Universidad de NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications 15 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-14483-x Sevilla, Seville, Spain. 42Fundació Docència i Recerca Mútua de Terrassa and Movement Disorders Unit, Department of Neurology, University Hospital Mutua de Terrassa, Terrassa, Barcelona, Spain. 43Hospital Universitario Central de Asturias, Oviedo, Spain. 44Hospital Universitario Parque Tecnologico de la Salud, Granada, Spain. 45Instituto de Investigación Sanitaria Biodonostia, San Sebastián, Spain. 46Hospital General de Segovia, Segovia, Spain. 47Memory Unit, Department of Neurology, IIB Sant Pau, Hospital de la Santa Creu i Sant Pau, Universitat Autònoma de Barcelona, Barcelona, Spain. 48Centro de Investigación Biomédica en Red en Enfermedades Neurodegenerativas (CIBERNED), Madrid, Spain. 49Hospital Universitario Marqués de Valdecilla-IDIVAL, Santander, Spain. 50University of Cantabria, Santander, Spain. 51Movement Disorders Unit, Department of Neurology, IIB Sant Pau, Hospital de la Santa Creu i Sant Pau, Universitat Autònoma de Barcelona, Barcelona, Spain. 52Centro de Investigacion Biomedica, Universidad de Granada, Granada, Spain. 53Hospital Universitario Virgen de las Nieves, Instituto de Investigación Biosanitaria de Granada, Granada, Spain. 54Hospital Clinic de Barcelona, Barcelona, Spain. 55Instituto de Investigación Sanitaria Fundación Jiménez Díaz, Madrid, Spain. 56Hospital Universitario Virgen de la Victoria, Malaga, Spain. 57Institut de Recerca Sant Joan de Déu, Barcelona, Spain. 58 Hospital Universitario Ramón y Cajal, Madrid, Spain. 59Department of Neurology, Instituto de Investigación Sanitaria La Fe, Hospital Universitario y Politécnico La Fe, Valencia, Spain. 60Hospital General de Segovia, Segovia, Spain. 61Department of Neurology, Hospital Universitario Fundación Alcorcón, Madrid, Spain. 62Department of Neurology, Medical University of Vienna, Vienna, Austria. 63Department of Neurology, Oslo University Hospital, Oslo, Norway. 64Department of Pathophysiology, University of Tartu, Tartu, Estonia. 65Department of Reproductive Biology, Estonian University of Life Sciences, Tartu, Estonia. 66Perron Institute for Neurological and Translational Science, Perth, WA, Australia. 67Department of Neurology and Neurosurgery, University of Tartu, Tartu, Estonia. 68Institute of Clinical Medicine, Department of Neurology, University of Oulu, Oulu, Finland. 69Department of Neurology and Medical Research Center, Oulu University Hospital, Oulu, Finland. 70University of Lago, Lagos State, Nigeria UK Brain Expression Consortium (UKBEC) Paola Forabosco71 & Robert Walker11 71 Istituto di Ricerca Genetica e Biomedica, Cittadella Universitaria di Cagliari, 09042 Monserrato, Sardinia, Italy 16 NATURE COMMUNICATIONS | (2020)11:1041 | https://doi.org/10.1038/s41467-020-14483-x | www.nature.com/naturecommunications