ARTICLE https://doi.org/10.1038/s41467-021-22262-5 OPEN Human-lineage-specific genomic elements are associated with neurodegenerative disease and APOE transcript usage Zhongbo Chen1,2,3, David Zhang1,2,3, Regina H. Reynolds 1,2,3, Emil K. Gustavsson 1,2,3, Sonia García-Ruiz1,2,3, Karishma D’Sa1,2,3, Aine Fairbrother-Browne1,2,3, Jana Vandrovcova1, International Parkinson’s Disease Genomics Consortium (IPDGC)*, John Hardy1,4,5,6,7, Henry Houlden 8, Sarah A. Gagliano Taliun 9,10, Juan Botía1,11 & Mina Ryten 1,2,3 ✉ 1234567890():,; Knowledge of genomic features specific to the human lineage may provide insights into brain- related diseases. We leverage high-depth whole genome sequencing data to generate a combined annotation identifying regions simultaneously depleted for genetic variation (constrained regions) and poorly conserved across primates. We propose that these con- strained, non-conserved regions (CNCRs) have been subject to human-specific purifying selection and are enriched for brain-specific elements. We find that CNCRs are depleted from protein-coding genes but enriched within lncRNAs. We demonstrate that per-SNP heritability of a range of brain-relevant phenotypes are enriched within CNCRs. We find that genes implicated in neurological diseases have high CNCR density, including APOE, highlighting an unannotated intron-3 retention event. Using human brain RNA-sequencing data, we show the intron-3-retaining transcript to be more abundant in Alzheimer’s disease with more severe tau and amyloid pathological burden. Thus, we demonstrate potential association of human- lineage-specific sequences in brain development and neurological disease. 1 Department of Neurodegenerative Disease, Queen Square Institute of Neurology, University College London (UCL), London, UK. 2 NIHR Great Ormond Street Hospital Biomedical Research Centre, University College London, London, UK. 3 Department of Genetics and Genomic Medicine, Great Ormond Street Institute of Child Health, University College London, London, UK. 4 Reta Lila Weston Institute, Queen Square Institute of Neurology, UCL, London, UK. 5 UK Dementia Research Institute, Queen Square Institute of Neurology, UCL, London, UK. 6 NIHR University College London Hospitals Biomedical Research Centre, London, UK. 7 Institute for Advanced Study, The Hong Kong University of Science and Technology, The Hong Kong University of Science and Technology, Hong Kong SAR, China. 8 Department of Neuromuscular Disease, Queen Square Institute of Neurology, UCL, London, UK. 9 Department of Medicine & Department of Neurosciences, Université de Montréal, Université de Montréal, Montréal, QC, Canada. 10 Montréal Heart Institute, Montréal, Québec, Canada. 11 Departamento de Ingeniería de la Información y las Comunicaciones, Universidad de Murcia, Murcia, Spain. *A list of authors and their affiliations appears at the end of the paper. ✉email: mina.ryten@ucl.ac.uk NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 H umans are perceived to be particularly vulnerable to 150 million single-nucleotide variants of which 54.7% had not neurodegenerative disorders relative to other primates on been reported in dbSNP16 or the most recent phase 3 of the 1000 both a pathological and phenotypic level1–5. This is Genomes Project17. The availability of this information has exemplified in Alzheimer’s disease, in which a similar phenotype already enabled more accurate identification of relatively con- is not seen in ageing non-human primates, nor are the char- strained regions of the genome, which has led to the development acteristic neurofibrillary tangles on pathological examination1,6. of the context-dependent tolerance score (CDTS)11. CDTS is Likewise, Parkinson’s disease does not naturally occur in non- derived from estimating how the observed genetic variation human primates, whose motor deficits do not respond to levo- compares to the propensity of a nucleotide to vary depending on dopa administration and a Lewy body pathological burden is not its surrounding context using the high-resolution profiles deter- present5,7. This has led to the hypothesis that the same evolu- mined from deep sequencing data11. Yet, this information has not tionary changes driving encephalisation which have steered the been combined directly with improved conservation data to development of characteristic human features may predispose to identify regions with evidence for human-specific selection. disorders that affect the brain2,5,6. In the case of Alzheimer’s In this study, we make full use of these resources to develop a disease, it is postulated that the accelerated evolution of intelli- novel, granular genomic annotation which efficiently captures gence, brain size and lifespan predisposes to selective advantages, information on intra-species constraint and inter-species con- which in later life have deleterious effects on cognition through servation simultaneously and identifies constrained, non- the very same pathways8. Therefore, identifying the genomic conserved regions (CNCRs). We use this annotation to test the changes unique to the human lineage may not only provide hypothesis that CNCRs are not only specific to the human line- insights into the evolution of human-lineage-specific phenotypic age, but given the encephalisation of humans, that CNCRs will be features but also into the pathophysiology underlying uniquely enriched within brain-specific functional and regulatory elements human diseases. as well as risk loci for neurological disease. We show that these Previous studies attempting to identify human-lineage-specific regions are enriched for SNP heritability for a range of neuro- variation and functional elements in the human genome have logical and psychiatric phenotypes. Furthermore, by calculating focused on genomic conservation as calculated by aligning and CNCR density within the boundaries of known genes, we develop comparing genomes across species. But, conservation measures a gene-based metric of human-specific constraint. This analysis alone do not fully identify regions with evidence of human- highlights APOE and leads to the identification of an intron-3 specific purifying selection. This is because a large part of the retaining transcript of APOE, the usage of which is correlated genome is evolving neutrally and sufficient phylogenetic distance with Alzheimer’s disease pathology and APOE-ε4 status. This is required to detect these changes9. Furthermore, alignment approach provides direct support for the role of human- methods do not reliably detect substitutions that preserve specific CNCRs in brain development and complex neurological function9. Conversely, some genes such as those implicated in phenotypes. immune system function may be subject to rapid evolutionary turnover even among closely related species9. For these reasons, analysing conservation alone has limited capacity to capture Results human-specific genomic elements9. Genomic regions with high constraint, but not conservation, The increasing availability of whole-genome sequencing were enriched for regulatory, non-coding genomic features. (WGS) has opened new opportunities to address this issue. Using CNC scores, which combine information from CDTS and intra-species whole-genome comparisons10,11, we are better able phastCons20, were used to capture evidence of disparity between to appreciate sequence differences between individuals of the constraint and conservation within a genomic region (Fig. 1). We same species, and identify genomic regions in humans containing investigated the relationship between CNC scores and known significantly fewer genetic variants than expected by chance, genomic features within the most constrained portion of the designated as constrained genomic regions. This form of analysis, genome (top 12.5%). This analysis demonstrated clear patterns of which is based on the assumption that most selection is negative enhancement and depletion for genomic elements across CNC or purifying (i.e. those that remove new deleterious mutations), scores, which significantly differed from similar analyses per- has been crucial for classification of exonic variation and attri- formed using constraint metrics alone11 (Fig. 2a). Among con- bution of pathogenicity12. However, many genomic regions strained genomic regions with the highest CNC scores (90 to 100 would be expected to be both constrained and conserved; such decile, signifying high constraint but low conservation), we saw a regions have been maintained by natural selection across species, depletion for coding elements of 27-fold relative to genomic including humans. This means that metrics reflecting constraint regions with the lowest CNC scores (chi-squared p < 2.2 × 10−16). alone cannot identify human-specific elements as the same This contrasts with the pattern using constraint metrics alone regions could also be conserved in other species. where the most constrained genomic regions are highly enriched This has led previous analyses to combine these metrics of for coding exons11. On the other hand, promoter, promoter- sequence constraint and conservation to identify genomic regions flanking and non-coding RNA features were overrepresented in with evidence for human-specific selection13,14. Ward and Kellis the highest compared to the lowest CNC deciles by 4.7- (chi- successfully applied this approach to demonstrate that a range of squared p < 2.2 × 10−16), 1.9- (chi-squared p < 2.2 × 10−16) and transcribed and regulatory non-conserved elements showed evi- 1.5-fold (chi-squared p < 2.2 × 10−16) respectively. Thus, genomic dence of lineage-specific purifying selection14. However, this regions with high CNC scores are enriched for regulatory, non- analysis was limited by the availability of WGS data and metrics coding genomic features. on human genetic variation were derived from the 1000 Genomes pilot data, which sequenced with only two to six times coverage15. Genes with the highest density of CNCRs are enriched for long Advances in sequencing technology have increased the feasibility non-coding RNA. Next, we applied a CNC score cut-off of ≥1 of deep sequencing of human populations leading to a much (signifying a twofold higher ranking in constraint than con- more detailed understanding of genetic variation between servation) to the 12.5% most constrained genomic regions to humans10. In fact, the recent sequencing of the genomes of define a set of genomic regions that were constrained, but not 10,545 human individuals at a coverage of 30–40 times identified conserved (termed CNCRs). We wanted to investigate whether 2 NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 ARTICLE Fig. 1 Workflow of study and schematic demonstration of annotation groups. The workflow depicts the processes involved in creation of the annotation with set parameters for each of the three groups of annotations generated and the processes involved in hypothesis-testing. CNC scores: constrained, non- conserved scores; CNCRs: constrained, non-conserved regions: CNCRs are defined as genomic regions that were first among the 12.5% most constrained, then with a CNC score of ≥1 (i.e. a twofold higher ranking in constraint than conservation). Constrained regions are defined as the regions within the 12.5% most constrained of the genome irrespective of conservation score. Non-conserved regions are defined as relatively non-conserved genomic regions with a conservation rank determined by the rank of the first quartile phastCons20 score at a CNC score of 1 (rank ≤ 25,623,592) (irrespective of constraint score). CDTS is the context-dependent tolerance score. Minus CDTS score is used as a lower score of CDTS corresponds to a more constrained region. CNCRs could be used to identify specific genes of interest. With Significant enrichment of heritability for neurologically rele- this in mind, we used CNCR density, the proportion of CNCRs vant phenotypes. Given the enrichment of regulatory features within a gene (defined in Supplementary Fig. 1), to identify gene within genomic regions with a high CNC score, we postulated sets which might be expected to contribute most to human- that such regions could also be enriched for disease risk. In order specific phenotypes. Consistent with the findings above, we found to study this, we investigated CNCRs for evidence of enriched that as the CNCR density threshold was increased to define heritability for a range of complex neurologically relevant phe- the gene sets of interest, there was a marked reduction in the notypes (Supplementary Table 4). After Bonferroni correction for proportion of protein-coding genes (β-coefficient between pro- multiple testing, we found that CNCRs exhibited significant portion and CNCR density = −1.061 and false discovery rate enrichment in heritability for intelligence test performance (FDR)-corrected p = 0.00162), and an increase in the proportion (coefficient p = 4.19 × 10−24); Parkinson’s disease (coefficient of long non-coding RNA (lncRNA, β-coefficient 0.385 and FDR- p = 4.65 × 10−5); major depressive disorder (coefficient p = correctedp = 0.0161) and microRNA-encoding genes (miRNA, 2.95 × 10−8) and schizophrenia (coefficient p = 5.26 × 10−19), but β-coefficient 0.394 and FDR-corrected p = 0.00116) (Fig. 2b). not for Alzheimer’s disease (Fig. 3). While a significant enrich- Interestingly, this relationship was not clearly observed when ment in heritability for intelligence test performance, major considering unprocessed snRNA and other RNAs (Fig. 2b). In depressive disorder and schizophrenia were also observed in the order to determine whether the relationship between CNCR constrained regions alone (and to a lesser extent, non-conserved density and gene biotype was driven by sequence constraint or regions), we noted that the regression coefficient for CNCRs was conservation, we also generated comparator gene lists based on at least twofold larger for the CNCR annotation compared to the constrained-only and non-conserved regions alone. Importantly, constrained annotation (Supplementary Table 4). Similarly, sig- lncRNA and protein-coding gene proportions do not follow the nificant enrichment in heritability for Parkinson’s disease was same directionality with increasing density when constraint or only observed in CNCRs. SNP heritability for Alzheimer’s disease non-conservation alone is considered (Fig. 2b). Thus, this analysis did not show significant enrichment although there was a trend highlighted the specific importance of lncRNAs as compared to for enrichment in terms of the regression coefficient and coeffi- other classes of non-coding RNAs in driving human-specific cient p value within CNCRs. Thus, by combining metrics for both patterns of gene expression. constraint and conservation in our annotation, we derived an NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 Fig. 2 Genomic territory and biotype proportions of constrained, non-conserved regions. Composition of the constrained genome, partitioned by constrained, non-conserved (CNC) scores (a) and proportion of biotypes of genes in our annotation (constrained, non-conserved regions: CNCRs) and in the comparator annotations (constrained regions and non-conserved regions) (b). The description for each genomic feature is shown in Supplementary Table 1. The barplot in a shows the genomic features for the 12.5% most constrained regions with CNC scores partitioned by decile, such that the highest decile (90–100) represents the most constrained and least conserved regions. Description of gene biotypes in b is taken from Ensembl42. The heatmap demonstrates the proportion of genes of a certain biotype within the three separate annotations within each genic CNCR density cut-off. CNCR density is defined as the proportion of CNCRs within a gene taking into account the gene size. Protein coding is defined by a gene that contains an open reading frame. The subclassified components of long non-coding RNA (lncRNA) found in the annotations are: Antisense—has transcripts that overlap the genomic span (i.e. exon or introns) of a protein-coding locus on the opposite strand; lincRNA (long interspersed ncRNA)—has transcripts that are long intergenic non-coding RNA locus with a length > 200 bp; non-coding RNA is further subclassified into miRNA (microRNA); siRNA (small interfering RNA); snRNA (small nuclear RNA) and miscellaneous RNA (includes snoRNA (small nucleolar RNA) and tRNA (transfer RNA)). Pseudogenes are similar to known proteins but contain a frameshift and/or stop codon(s) which disrupts the open reading frame. These can be classified into processed pseudogene—a pseudogene that lacks introns and is thought to arise from reverse transcription of mRNA followed by reinsertion of DNA into the genome and unprocessed pseudogene—a pseudogene that can contain introns since produced by gene duplication. independent annotation that shows a higher level of enrichment evolutionary pressures, the analysis of gene sets defined in heritability for neurologically related phenotypes than anno- on the basis of non-conservation alone demonstrated the tations based on constraint or conservation alone. significant enrichment of immune and skin-related GO terms (GO:0002250 adaptive immune response, p = 4.02 × 10−10 and GO:0043588 skin development, p = 2.33 × 10−4). Taken together, The proportion of enriched gene sets with neurologically these results demonstrate that using CNCR density, genes related GO terms increases in genes with the highest density of important in nervous system development and implicated in CNCRs. To investigate these findings further, we defined gene neurological disease can be identified. sets based on their CNCR density (the proportion of CNCRs within a gene) and analysed their GO term enrichment. We assessed gene sets defined across a range of CNCR densities (>0.0 CNCR annotation highlights an intron-3 retaining transcript to ≥0.5 at 0.1 increments). We found that the proportion of of APOE. Next, we investigated the distribution of CNCR density neurologically associated GO terms with significant enrichments across Mendelian genes associated with a neurological pheno- (g:SCS-corrected p < 0.05) increased among gene sets with types (as defined within Online Mendelian Inheritance in Man increasing CNCR gene densities (Supplementary Fig. 2). Impor- (OMIM)18) and genes implicated in complex brain-relevant tantly, a similar analysis of gene sets defined by constraint alone phenotypes (as defined within Systematic Target OPportunity or non-conservation alone did not contain any neurologically assessment by Genetic Association Predictions (STOPGAP)19). enriched GO terms (Fig. 4). We identified the gene set with the We noted that the median CNCR density was significantly highest proportion of nervous system-related terms at a CNCR higher in OMIM genes with a neurological phenotype compared genic density of 0.3 (Supplementary Fig. 2). The only GO terms to all other genes (median CNCR density of neurological OMIM specific to a tissue process were related to the nervous system genes = 0.0924, IQR = 0.0567 − 0.143; median CNCR density of (Fig. 4, Supplementary Table 5) and spanned terms such as all other genes = 0.083, IQR = 0.043 − 0.153; Wilcoxon rank sum neuronal development (GO:0048663, corrected p = 5.46 × 10−7) test p = 1.8 × 10−6). While genes associated with complex brain- and spinal cord differentiation (GO:0021515, corrected p = relevant phenotypes did not have a significantly higher CNCR 3.64 × 10−7). The remaining significantly enriched GO terms density when compared to all other genes, we still identified 31 related to ubiquitous processes including protein targeting genes with a CNCR density of greater than 0.2 and 7 genes with a (GO:0045047, p = 9.93 × 10−4) and DNA binding (GO:0043565, CNCR density of greater than 0.3 (APOE, PHOX2B, SSTR1, p = 4.81 × 10−4). Of note, analysis of gene sets defined on the HCFC1, HAPLN4, CENPM and IQCF5). Of these genes, APOE basis of constraint alone revealed no enrichment of neurologically had the highest CNCR density with a value of 0.552. associated terms, but instead significant enrichment of vascular Given the high CNCR density of APOE, its importance as a system-related GO terms (GO:0048514 blood vessel morpho- disease locus for Alzheimer’s disease and other neurodegenerative genesis, corrected p = 3.96 × 10−37 and GO:0072358 cardiovas- diseases20 and the long-standing interest in its lineage cular system development, p = 8.53 × 10−36). As might be specificity8,21 (specifically the differences in the ɛ4 allele between expected based on the rapid and potentially divergent humans and non-human primates1), we chose to focus on this 4 NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 ARTICLE Fig. 3 Stratified-linkage disequilibrium score regression (s-LDSC) analysis across five traits comparing constrained, non-conserved regions (CNCRs) with its constituent constrained and non-conserved annotations. a The regression coefficient. b The regression coefficient −log10(p value) with the dotted line showing the Bonferroni-corrected p value of 0.00333 for 15 conditions. GWASs were as follows: Intelligence2019: intelligence test performance GWAS, AD2018: Alzheimer’s disease GWAS, PD2019.ex23&Me: Parkinson’s disease GWAS without 23&Me data, MDD2018: major depressive disorders GWAS and SCZ2018: schizophrenia GWAS (Supplementary Table 2). Fig. 4 Summarised enriched gene sets for terms specific for neurological gene sets, other non-neurological-specific tissues and non-tissue-specific as defined by Gene Ontology (GO). Plot comparing annotation of interest (CNCRs) and comparator annotations which only use constraint or non-conserved metrics. Frequency, derived from REViGO49, represents the percentage of human proteins in UniProt which were annotated with a GO term, i.e. a higher frequency denotes a more general term. NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 Fig. 5 Annotation with constrained, non-conserved regions (CNCRs) is highly granular and shows APOE to have a high density of CNCRs throughout its length especially in association with an intron-3 retention event in the human hippocampus. The first track represents the genomic location of APOE within chromosome 19. The second track shows the known transcripts, currently within annotation in Ensembl v.92. The mean coverage (MC) (log10 scale) in the hippocampus shown here is greater than zero (denoted by the grey shaded area) across intron-3 highlighting an intron-3 retention event (mean coverage data derived from GTEx v.7). In the fourth track, CNC scores above the black dashed line and shaded in red fulfil criteria for a constrained, non- conserved region (CNCR) are shown. The intron-3 retention event has the highest CNCR density among all intronic regions of APOE. The fifth track labelled “CNCR” depicts regions fulfilling criteria for CNCR. PC20 represents the mean phastCons20 score. The black dashed line within this track represents a mean phastCons20 score of 0. CDTS represents the context-dependent tolerance score as a measure of constraint with the black dashed line showing a CDTS of 0. Within the CDTS track, the blue dotted line represents a region with no CDTS annotation. gene to further validate our annotation. We tested whether observed for the three human transcriptomes but do not occur in intragenic analysis of APOE could identify specific regions or abundance for the three chimpanzee samples (Supplementary transcripts of interest. We compared CNCR density, constraint Fig. 5b, c). We also noted that there was a trend for lower and conservation scores across the length of the gene showing expression of intron-3 within chimpanzees compared to humans that CNCRs provide a highly granular annotation (Fig. 5). Using (Supplementary Fig. 6). However, this analysis was limited by the this approach, we identified the region with the highest CNCR inherently small sample sizes and so comparing the coverage of density in APOE to be within intron-3 (Supplementary Fig. 3), intron-3 normalised for the total coverage of APOE within the coinciding with the annotated region’s boundaries. Furthermore, samples did not identify statistically significant differences in the the intron-3 region had a higher coverage compared to introns expression of intron-3 across species. Nonetheless, these data 1 and 2 based on the mean coverage provided by Genotype- would suggest that intron-3 retaining transcripts are more Tissue Expression Consortium (GTEx) hippocampal tissue commonly expressed in humans, and likely to be largely absent indicating that this is likely to represent an intron retention in chimpanzees. event (Supplementary Fig. 3E). This coverage was calculated as In order to obtain further insights into the biological the mean across all GTEx samples normalized to a target library significance of the intron-3 retaining APOE transcript, we size of 40 million 100 base pair (bp) reads (mean coverage seen in leveraged publicly available RNA-sequencing data covering 11 Fig. 5)22. Thus, in conjunction with the highest intragenic CNCR regions of the human central nervous system provided by the density localised to intron-3, these coverage data provided further GTEx v.725. Using an annotation-independent approach to justification for our analysis of the intron-3 retention event. identify genomic regions producing stable transcripts26,27, we Although no intron-3 retaining transcript is currently anno- identified a region of significant expression encompassing intron- tated in Refseq and Ensembl, an intron-3 retention event has 3 of APOE and the flanking coding exons in all brain tissues previously been reported and implicated in the regulation of (Fig. 6a). These data not only support the existence of an intron-3 APOE expression21,23,24. To validate the existence of this retaining APOE transcript that is not entirely attributable to pre- transcript, we performed Sanger sequencing of poly-A-selected mRNA transcripts or driven by background noise in sequencing RNA derived from human hippocampal tissue. This demon- but also provide a means of estimating its usage across the human strated the existence of a transcript containing the full-length brain. intron-3 sequence was flanked by both exon 3 and exon 4 Thus, in order to compare usage of this transcript across (Supplementary Fig. 4a). Using a non-RT control, we showed that different CNS regions, we calculated the ratio of normalised this could not be explained by genomic DNA (gDNA) intron-3 expression (a measure of intron-3 retaining transcripts) contamination (Supplementary Fig. 4b). to the normalised expression of exon 3/exon 4 spanning reads We noted that transcripts retaining intron-3 of APOE are (a measure of transcripts splicing out intron-3). We see that there unannotated by Ensembl for humans and chimpanzees (Supple- is evidence of the usage of the intron-3 retaining APOE transcript mentary Fig. 5a). However, reads aligning to intron-3 were in all central nervous system regions from GTEx data (Fig. 6a). 6 NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 ARTICLE Fig. 6 Quantification of APOE intron-3 retaining transcript usage. Quantification of intron retention usage by its normalised coverage to junction ratio across brain tissues within GTEx (a). Normalised coverage to junction ratio of the APOE intron-3 retention event in bulk RNA-sequencing data of post- mortem dorsolateral prefrontal cortex tissue samples from 634 individuals recruited within ROSMAP studies across Braak and Braak staging (b) and APOE ɛ4 allele status (c). In a, red dashed horizontal line presents the median normalised intron retention coverage to junction ratio within central nervous system tissues in GTEx. Number of samples within each of the tissue groups was as follows: amygdala—72; anterior cingulate cortex—84; caudate—117, cerebellar hemisphere—105; frontal cortex—108; hippocampus—94; hypothalamus—96; nucleus accumbens—113; putamen—97; spinal cord—71; substantia nigra—63. The Kruskal–Wallis p value show results from comparison of the differences in the normalised intron retention coverage to junction ratio between the different brain regions with pairwise regions comparisons shown in Supplementary Table 6. In b and c, the blue line represents the linear regression fit with the grey shaded area representing ± 95% confidence interval. Braak and Braak staging is a measure of severity of neurofibrillary tangle based on location. To improve the power of the study, we merged Braak and Braak stages I and II to “Braak mild stage”, Braak and Braak stages III and IV to “Braak moderate” and Braak and Braak stages V and VI to indicate “Braak severe” stage. For number of APOE ɛ4 alleles, a heterozygous state is represented by “1” and homozygous state by “2”. However, there are also significant differences among brain Memory and Aging Project (ROSMAP) studies to quantify the regions (Kruskal–Wallis p < 2.2 × 10−16) (Supplementary Table 6: usage of the intron-3 retaining transcript of APOE in post- pairwise Wilcoxon rank sum test p values) with the usage of the mortem dorsolateral prefrontal cortex brain tissue derived from intron-3 retaining event being highest in the spinal cord, individuals with Alzheimer’s disease (n = 222) and mild cognitive substantia nigra and hippocampus (Fig. 6a). impairment (MCI) (n = 158) compared to control individuals In summary, we confirmed the existence of an unannotated (defined as the final clinical diagnosis blinded to pathological human-specific non-coding transcript of APOE and identified findings, n = 202). Prior to our analyses, we assessed the impact differential usage of this transcript across the human brain. In this of batch effects within this dataset. After finding that our analyses way, we demonstrated the utility of combining CNC scores with were robust to the removal of an outlying batch (batch 7, transcriptomic data, which we have made easier through the Supplementary Fig. 7), we incorporated all batches into the visualisation platform vizER (https://snca.atica.um.es/browser/ analyses. Using this approach, we found that the proportion of app/vizER). Furthermore, this direct visualisation allows identi- the intron-3-retaining transcript was higher (p < 2.2 × 10−16) in fication of isolated intragenic regions of functional importance in dorsolateral prefrontal cortex tissue from individuals with genes with highly variable CNCR density. clinically diagnosed Alzheimer’s disease and MCI patients versus control participants. Partitioning this further on the basis of Usage of the intron-3 retaining transcript of APOE correlates pathology, we see an increase in intron-3 retaining transcript with Alzheimer’s disease pathology and APOE genotype. We usage with more severe Braak and Braak pathology for noted that among the brain tissues with the highest usage of the neurofibrillary tangles (adjusted r2 0.678, p < 2.2 × 10−16) intron-3 retaining transcript of APOE are those that show (Fig. 6b). Consistent with these findings, we also found a selective vulnerability for neurodegeneration, namely the hippo- significant increase in transcript usage with higher amyloid campus in the context of Alzheimer’s disease, the substantia nigra plaque pathology as defined using CERAD staging (adjusted r2 in the context of Parkinson’s disease and the spinal cord in the 0.673, p < 2.2 × 10−16). Finally, we investigated the relationship context of amyotrophic lateral sclerosis (pairwise comparisons between presence of the ε4 allele in APOE and usage of the between brain regions shown in Supplementary Table 6). intron-3 retaining transcript. We found a significant positive Given that APOE is one of the most important genetic risk correlation between ε4 allele load and the proportion of intron-3 factors for Alzheimer’s disease, we leveraged publicly available retaining transcript (adjusted r2 0.673, p < 2.2 × 10−16) RNA-sequencing data from the Religious Orders Study and (Fig. 6c). This association remained significant after partitioning NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 APOE-ε4 status by disease and accounting for tau and amyloid of mammalian organ development32. We also find enrichment of burden, showing that this association is likely to be independent spinal cord-associated genes that may relate to the uniquely of disease state. human monosynaptic corticomotoneuronal pathways implicated Taken together, these findings could suggest that usage of the in human-specific dexterity and digital motor control33,34, the intron-3 retaining transcript may be regulated by APOE-ε4 status disruption of which may lead to amyotrophic lateral sclerosis35. and may be involved in mediating the effect of APOE genotype, We noted that APOE was among the genes with the highest supporting a role for the presence of this lncRNA in disease risk CNCR density across the genome and carried the highest CNCR and progression, although it is also feasible that Alzheimer’s density of all genes implicated in complex brain-relevant phe- disease pathology could drive intron-3 retention notypes (defined within the STOPGAP database19). Given that genetic variation within this gene and specifically APOE-ε4 status is not only the principal genetic risk factor for Alzheimer’s Discussion disease36 but also associated with risk for other neurodegenerative The core aim of this study was to test the hypothesis that cap- disorders, stroke and reduced lifespan20, this finding provides turing human-lineage-specific regions of the genome could pro- evidence for the value of CNCR annotation. We thus further vide insights into neurological phenotypes and diseases in studied APOE to validate our annotation. Within APOE, the humans. We generated and used an annotation based on existing CNCR annotation highlighted an intron-3 retention event of high knowledge of sequence conservation and sequence constraint coverage and CNCR density, not currently within annotation but within humans, which we termed CNCRs. We used this anno- which has been previously reported to be associated with neu- tation to prioritise genomic regions, genes and transcripts based ronal regulation of APOE, with splicing out of the intron-3- on a high density of human-lineage-specific sequence as deter- containing mRNA following neuronal injury in neuronal cell lines mined by our CNCR annotation. We demonstrated the utility of and human APOE knock-in mouse models21,23,24. Using Sanger this approach by showing: the genomic regions we identified are sequencing of cDNA derived from control human hippocampal enriched for SNP heritability for intelligence test performance tissue, we confirm the presence of an intron-3 retaining APOE and brain-related disorders; the genes we identified are enriched transcript. We estimated the usage of the transcript from short- for neurologically relevant gene ontology terms and genes causing read RNA-sequencing data and found variable levels across dif- neurogenetic disorders and the existence of an intron-3 retaining ferent brain tissues within GTEx25 with the highest usage in the transcript of APOE, the usage of which is correlated with Alz- spinal cord, substantia nigra and hippocampus, reflecting central heimer’s disease pathology and APOE-ε4 status. nervous system regions most susceptible to selective vulnerability A major finding of this study is that CNCRs are enriched for in disease. Using human dorsolateral prefrontal cortex RNA- regulatory, non-coding genomic regions. This is consistent with sequencing data, we found that the intron retention event was analyses performed by Ward and Kellis14, and highlights the significantly more abundant in patients with Alzheimer’s disease potential functional importance of non-conserved and thus evo- than controls and in those with more severe Braak and Braak lutionarily recent non-coding regions subject to constraint. Fur- pathology and amyloid burden as characterised by CERAD thermore, these findings suggest that CNCRs could provide a pathology. Furthermore, we saw a dosage-dependent increase in means of prioritising and potentially aiding the assessment of the intron retention event with the APOE-ε4 allele that was non-coding variants, an area of significant interest, given that independent of disease status. Although our findings do not 88% of GWAS-derived disease-associated variants reside in non- elucidate the function of the intron-3 retention event, they are coding regions of the genome28. We found evidence to support consistent with previous studies that have shown general this view through heritability analyses for intelligence test per- increases in intron retention events as a feature of Alzheimer’s formance, Parkinson’s disease, major depressive disorder and disease and ageing with implications for post-transcriptional schizophrenia with SNP heritability not only enriched within regulation37. We propose that this novel transcript may be a CNCRs, but to a greater extent than would be expected using means of regulating APOE in a disease state or could itself be either conservation or constraint annotations alone. Considering driven by Alzheimer’s disease pathology. heritability for intelligence test performance, this phenotype is Given that we use existing measures of constraint and con- already known to also be enriched within annotations of brain- servation to identify CNCRs, this analysis is fundamentally lim- specific tissue expression and among several regulatory biological ited by the quality of these data. While the constraint metrics we gene sets29, including neurogenesis, central nervous system neu- used were derived from high-depth sequencing, this is still ron differentiation and regulation of synapse structure or restricted given the relatively high number of private genetic activity28. These findings support our hypothesis that CNCRs variants we each carry. In addition, analysis was limited to identify genomic regions of functional importance with relevance the high-confidence regions covering ~84% of the genome, to human brain phenotypes. amounting to 12.2% of all genes that remained unannotated with Our analyses of CNCR density within genes are consistent with CDTS metrics11. Thus, on balance, it is difficult to predict the these findings, highlighting both non-coding genes and those impact of these missing data on our findings. Similarly, our study implicated in neurologically relevant processes and diseases. of the relationship between CNCRs and known genomic features Interestingly, CNCR annotation specifically highlighted lncRNAs is limited by the annotation quality in existing databases. We have as opposed to other non-coding RNAs. In particular, we observed endeavoured to overcome some of these problems by creating a a proportional increase in lncRNA enrichment with higher genic more detailed annotation combining both GENCODE and CNCR density, which could not be replicated using measures of Ensembl data as used by di Iulio et al. in their work generating sequence constraint or conservation alone. This observation is in CDTS11. The SNP heritability estimates using stratified-linkage keeping with previous studies that have shown most lncRNAs are disequilibrium score regression (LDSC) analysis are limited by tissue-specific with the highest proportion being specific to the quality of linkage disequilibrium (LD) information under- brain30 and highly relevant to neurodegenerative diseases31. pinning the heritability calculations38 and the sample size of Similarly, the enrichment for nervous system-related pathways the GWAS. within CNCRs, which is representative of recent purifying Despite these limitations, we have been able to demonstrate the selection, is in keeping with the lowest proportion of positively utility of CNCRs specifically in the identification of functionally selected genes being present in brain tissues from previous studies important non-coding regions of the genome, genes and 8 NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 ARTICLE transcripts. We find that CNCRs across all forms of analyses (Fig. 1) which fall either into constrained-only or non-conserved only annotations highlight the significance of human-lineage-specific sequences in as defined respectively by: (i) CNCRs annotation: SNPs falling into CNCRs; (ii) constrained annotation: SNPs located within the 12.5% most constrained regions of the central nervous system and in the context of neurological the genome irrespective of conservation score and (iii) non-conserved annotation: phenotypes and diseases. We release our annotation of CNC SNPs located within relatively non-conserved genomic regions with a conservation scores and CNCRs via the online platform vizER (https://snca. rank determined by the rank of the first quartile phastCons20 score at a CNC score atica.um.es/browser/app/vizER) to allow CNCRs to be viewed at a of 1 (rank ≤ 25,623,592) (irrespective of constraint score). We provided Bonferroni- granular level. Thus, the CNCR annotation we generate has the correctedp values, which account for the number of annotation categories and GWASs tested (total of 15 conditions). potential to provide additional disease insights beyond those explored within this study and as we anticipate the release of increasing quantities of WGS data in humans will only improve Generation of a gene-based metric for CNCRs and gene set enrichment analysis. To generate a metric of human-specific constraint, which could be in quality and value. applied to a gene rather than a 10 bp region, we calculated the density of CNCRs within each gene, the length of which was defined by the transcription start and Methods stop sites for that gene (GRCh38.v97). The CNCR density was defined as the Generation of an annotation for the identification of CNCRs. We generated a proportion of the length of a gene containing CNCRs (Supplementary Fig. 1d). In combined annotation to capture information on intra-species constraint and inter- this way, we were able to normalise for the effect of gene size on our metric. species conservation simultaneously, using CDTS together with phastCons20 scores Therefore, the higher the gene density, the larger the proportion of the total length (Fig. 1). The previously validated map of sequence constraint (http://www.hli- of the gene was covered by CNCRs. opendata.com/noncoding) generated using 7794 whole-genome sequences11 was In order to compare the relationship between the change in CNCR density and used to assign a single CDTS score to each non-overlapping 10 bp region throughout the proportion of a genic biotype (defined by Ensembl v.92), we used linear the genome (build GRCh38, 248,925,226 bins). The phastCons20 score, which cal- regression and applied FDR-correctedp values in R v.3.6.1. culates the likelihood ratio of negative selection based on the total number of sub- We used g:ProfileR (R Package)48 for gene set enrichment analysis. We used the stitutions during evolution of an element between species39, was used as a measure of three sets of tested annotations incorporating genes that fell into CNCRs, inter-species conservation (http://hgdownload.cse.ucsc.edu/goldenPath/hg38/ constrained regions and non-conserved regions in the gene set enrichment analysis phastCons20way/)39. PhastCons20 was used as it compares the human genome to as previously described for LDSC annotation and as defined in Fig. 1. The the genomes of less divergent species (16 other primates and three mammals). For background gene list in all analyses comprised 49,644 genes from all regions of the each 10 bp bin, we assigned the corresponding mean phastCons20 score. Bins genome with a CDTS and phastCons20 annotation. The correction method was set without a conservation score due to insufficient species in the alignment were not to g:SCS to account for multiple testing48. We used REViGO49 to summarise the considered (0.218% of the genome), nor did we consider bins without a CDTS score significant GO terms, and to derive the term frequency, which is a measure of GO (16% of the genome, equating to 12.2% of all genes). We found that 10.9% of the term specificity. unannotated regions of the genome were within the ENCODE list of problematic To further characterise CNCR density within genes associated with disease, we regions (https://github.com/Boyle-Lab/)40 with the remainder accounted for by first studied phenotype relationships of all Mendelian genes within the OMIM incomplete sequencing from the 10,000 Genomes Project10. We recognise that this is catalogue (http://api.omim.org)18. We compared the CNCR density of all a limitation of this study and of the previously reported analysis11. For the remaining neurologically relevant OMIM genes to all genes within CNCR annotation. 248,381,744 bins, we ranked both CDTS and mean phastCons20 scores across the Secondly, in order to investigate the CNCR density within genes associated with whole genome such that the highest ranks represented the most constrained and complex disorders, we used the STOPGAP database, a catalogue of human genetic conserved regions, respectively. We calculated the log2 ratio of the rank of constraint associations mapped to effector gene candidates derived from 4684 GWASs19. We to the rank of conservation for each 10 bp bin (termed constrained, non-conserved selected for genes associated with SNPs that surpassed a genome-wide significant score, CNC score). This resulted in scores with a distribution centred at 0 signifying p value of 5 × 10−8 and which fulfilled medical subject heading for associated no fold change between the ranks of the two metrics (Supplementary Fig. 1). Finally, neurological/behavioural diseases. We used these sets to identify potential genes of we defined CNCRs as genomic regions that were first among the 12.5% most interest associated with brain-related disorders which carry a high CNCR density. constrained, then with a CNC score of ≥1 (i.e. a twofold higher ranking in constraint than conservation). We used this definition for CNCRs throughout this study to Sequencing of APOE transcripts in human brain. Focussing on a region with capture regions that were among the most constrained, but less conserved genome. high CNCR density identified within APOE from the preceding analyses, we used Sanger sequencing of cDNA reverse transcribed from pooled human hippocampus Investigating the relationship between CNCRs and existing annotation. To poly-A-selected RNA (Takara/Clontech 636165) to support the presence of the investigate the relationship between CNC scores for genomic regions and genomic intron-3 retention event identified within APOE (GRCh38: chr19:44907952- features, we calculated the distribution of CNC scores across genomic features 44908531). For the reverse transcription, we used 500 ng of input RNA, with defined by GENCODE v.5341 and Ensembl v.9242. We restricted our analysis to the 10 mM dNTPs (NEB N0447S), VN primers and strand-switching primers (Oxford 12.5% most constrained regions only (31,115,616 ten bp bins) and segregated these Nanopore Technologies SQK-DCS109), 40 units of RNaseOUT inhibitor (Life regions into equally sized deciles ranked on the basis of CNC scores such that the Technologies 10777019) and 200 units of Maxima H Minus reverse transcriptase highest decile (90–100 decile) represented a high CNC score containing the most with 5X reverse transcription buffer (Thermo Fisher EP0751). PCR amplification of constrained and least conserved sequences. Each 10 bp region was then assigned a the cDNA was performed using primer pairs designed to span across intron-3 and single overlapping genomic feature. To avoid conflicts arising from overlapping exon 4 (P2-4) and intron-3 alone (P5) of APOE (ENST00000252486.9) (Supple- GENCODE and Ensembl definitions, we preferentially assigned a single genomic mentary Table 3). PCR was performed using Taq DNA polymerase with Q- feature to a given region by prioritising features as used by di Iulio et al.11 solution (Qiagen) and enzymatic clean-up of PCR products was performed using (described in Supplementary Table 1). In order to compare the enrichment of Exonuclease I (Thermo Fisher Scientific) and FastAP thermosensitive alkaline existing annotations within the proportions of the different genic regions, we used phosphatase (Thermo Fisher Scientific). Sanger sequencing was performed using chi-squared test with Yate’s continuity correction, implemented in R v.3.6.1. the BigDye terminator kit (Applied Biosystems) and sequence reactions were run on ABI PRISM 3730xl sequencing apparatus (Applied Biosystems). Electro- pherograms were viewed and sequences were exported using Sequencher 5.4.6 Enrichment of common-SNP heritability in brain-related phenotypes for (Gene Codes). Sequences were aligned against the human genome (hg38) using CNCRs. Stratified-LDSC was used to assess the enrichment of common-SNP BLAT and visually inspected for confirmation of validation. heritability for a range of complex diseases and traits within our annotation38,43. To reduce the risk of gDNA contamination, the human hippocampus poly-A- Stratified-LDSC makes use of the increased likelihood of a causal relationship in a selected RNA (Takara/Clontech 636165) had undergone selection by two rounds of block of SNPs in LD to correct for confounding biases that include cryptic relat- oligo(dT)-cellulose columns. Furthermore, the Maxima H Minus reverse transcriptase edness and population stratification in a polygenic trait43. Using established pro- buffer contains a double-strand-specific DNase to specifically remove gDNA. Lastly, tocols (https://github.com/bulik/ldsc/wiki), we tested whether SNPs located within we used poly-A selected RNA sample as a no-reverse transcriptase control in our annotation contributed significantly to SNP heritability after controlling for a comparison with cDNA to show that there is no contamination with gDNA. A total range of other annotations described within the baseline mode (v.1.2). This analysis of 100 ng poly-A-selected RNA and 100 ng cDNA was used for each reaction in a generates a coefficient z-score, from which we calculated a one-tailed coefficient total volume of 10 µl using the same PCR conditions. PCR amplification of cDNA and p value. Stratified-LDSC regression analyses were also run to incorporate back- RNA was performed using primer pair P2 using Taq DNA polymerase (Qiagen ground SNPs defined as all SNPs in the genome that include a CDTS and phast- protocol) with 30 s of denaturation and 30 s of annealing at 57 °C. Cons20 annotation, to avoid overestimation of the contribution to SNP heritability. We assessed the annotation for SNP heritability enrichment in complex brain- related disorders and phenotypes of intelligence test performance29, Alzheimer’s Quantifying differences in intron-3 retention between different species. In disease44, Parkinson’s disease (excluding 23&Me participants)45, schizophrenia46 order to investigate species differences in the human APOE intron-3 retention and major depressive disorder (excluding 23&Me participants)47 (Supplementary event identified with chimpanzees, we leveraged existing bulk RNA-sequencing Table 2). We considered SNPs within CNCRs and its two constituent groups data derived from chimpanzee and human hippocampus reported in Khrameeva NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 et al.50. We used data from the hippocampus and downloaded FASTQ files (NCBI 2. O’Bleness, M., Searles, V. B., Varki, A., Gagneux, P. & Sikela, J. M. Evolution Gene Expression Omnibus; https://www.ncbi.nlm.nih.gov/geo/, accession number of genetic and genomic features unique to the human lineage. Nat. Rev. Genet. GSE127898) from the three available human samples (SAMN11165674, 13, 853–866 (2012). SAMN11165673, SAMN11165737) and three available chimpanzee samples 3. Xu, K., Schadt, E. E., Pollard, K. S., Roussos, P. & Dudley, J. T. Genomic and (SAMN11166008, SAMN11165613, SAMN11165949). All bulk RNA-sequencing Network patterns of schizophrenia genetic variation in human evolutionary FASTQ files were aligned using STAR51 with a reference index generated for the accelerated regions. Mol. Biol. Evol. 32, 1148–1160 (2015). relevant species. Given the small number of samples, direct visualisation of the 4. Cookson, M. R. Evolution of neurodegeneration. Curr. Biol. 22, R753–R761 aligned BAM files in Integrative Genome Viewer52 across APOE was carried out. (2012). Furthermore, for all samples, the total coverage of intron-3 was normalised for the 5. Diederich, N. J., James Surmeier, D., Uchihara, T., Grillner, S. & Goetz, C. G. length of intron-3 in each species and also for the total APOE coverage accounting Parkinson’s disease: is it a consequence of human brain evolution? Mov. for sequencing depth and APOE expression differences to allow cross-species and Disord. 34, 453–459 (2019). cross-sample comparisons. The co-ordinates for the regions were taken from 6. Gearing, M., Rebeck, G. W., Hyman, B. T., Tigges, J. & Mirra, S. S. Ensembl for human GRCh38 and chimpanzee Pan_tro_3.0. Neuropathology and apolipoprotein E profile of aged chimpanzees: implications for Alzheimer disease. Proc. Natl Acad. Sci. USA 91, 9382–9386 Analysis of public RNA-sequencing data. We used publicly available short-read (1994). RNA-sequencing data from human brain post-mortem samples provided by GTEx 7. Collier, T. J., Kanaan, N. M. & Kordower, J. H. Ageing as a primary risk factor v.7.125 and the ROSMAP Study53 and to quantify the intron-3 retention event in for Parkinson’s disease: evidence from studies of non-human primates. Nat. APOE highlighted by our analysis. For GTEx data, we used pre-aligned files Rev. Neurosci. 12, 359–366 (2011). available from recount2 (https://jhubiostatistics.shinyapps.io/recount/)54. Both 8. Raichlen, D. A. & Alexander, G. E. Exercise, APOE genotype, and the studies within ROSMAP are longitudinal clinicopathological cohort studies of evolution of the human lifespan. Trends Neurosci. 37, 247–255 (2014). aging and/or Alzheimer’s disease. We downloaded BAM files for ROSMAP bulk 9. Kellis, M. et al. Defining functional DNA elements in the human genome. RNA-sequencing data from the Synapse repository (https://www.synapse.org/#! Proc. Natl Acad. Sci. USA 111, 6131–6138 (2014). Synapse:syn4164376) for analysis. To quantify the intron-3 retention event, we 10. Telenti, A. et al. Deep sequencing of 10,000 human genomes. Proc. Natl Acad. calculated the coverage of intron-3 expression normalised for the coverage across Sci. USA 113, 11901–11906 (2016). the entire APOE gene, as defined by the transcription start and end sites. To 11. di Iulio, J. et al. The human noncoding genome defined by genetic diversity. quantify splicing of intron-3, we calculated the number of exon 3 to exon 4 junc- Nat. Genet. 50, 333–337 (2018). tion reads (defined as reads mapping with a gapped alignment), normalised for all 12. Lek, M. et al. Analysis of protein-coding genetic variation in 60,706 humans. APOE junction reads detected and currently within annotation. We used a ratio of Nature 536, 285–291 (2016). the normalised coverage to normalised junction count over intron-3 as an estimate 13. Schrider, D. R. & Kern, A. D. Inferring selective constraint from population of the proportional use of the intron-3-retaining transcript, such that a high ratio is genomic data suggests recent regulatory turnover in the human brain. Genome associated with a higher usage of intron retention within both GTEx and ROSMAP Biol. Evol. 7, 3511–3528 (2015). data. Normalisation of the intron-3 event for APOE gene expression (directly 14. Ward, L. D. & Kellis, M. Evidence of abundant purifying selection in humans proportional to the canonical transcripts) was used to show independent effects of for recently acquired regulatory functions. Science 337, 1675–1678 (2012). the intron-3 event from the canonical transcripts. Comparisons between the two 15. The Genomes Project C. et al. A map of human genome variation from groups were performed by comparing the mean values of this normalised measure using Wilcoxon rank sum test, taking two-tailedp values < 0.05 to be significant. population-scale sequencing. Nature 467, 1061 (2010). Based on existing ROSMAP results55 and principal component analysis of frag- 16. Sherry, S. T. et al. dbSNP: the NCBI database of genetic variation. Nucleic ments per kilobase million data, we incorporated covariates to account for the Acids Res. 29, 308–311 (2001). effect of batch, RNA integrity number, post-mortem interval, study index, ethni- 17. Sudmant, P. H. et al. An integrated map of structural variation in 2,504 city, age at death and sex on estimates of intron-3-retaining transcript usage. Using human genomes. Nature 526, 75 (2015). the resulting linear regression model, we compared the intron-3 retention nor- 18. Hamosh, A., Scott, A. F., Amberger, J., Valle, D. & McKusick, V. A. Online malised coverage to junction ratio across clinical disease states, pathological states Mendelian Inheritance in Man (OMIM). Hum. Mutat. 15, 57–61 (2000). and APOE-ɛ4 status in 634 post-mortem brain samples. 19. 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). Reporting summary. Further information on research design is available in the Nature 20. Belloy, M. E., Napolioni, V. & Greicius, M. D. A quarter century of APOE and Research Reporting Summary linked to this article. Alzheimer’s disease: progress to date and the path forward. Neuron 101, 820–838 (2019). Data availability 21. Mahley Robert, W., Huang, Y. & Apolipoprotein, E. Sets the stage: response to We release our annotation of CNC score as an interactive visualisable track via online injury triggers neuropathology. Neuron 76, 871–885 (2012). platform vizER: (https://snca.atica.um.es/browser/app/vizER) and provide a publicly 22. Zhang, D. et al. Incomplete annotation has a disproportionate impact on our downloadable table of CNCR density for genes within our annotation (under the understanding of Mendelian and complex neurogenetic disorders. Sci. Adv. 6, “Download” Tab). eaay8299 (2020). Publicly available datasets used are: 23. Xu, Q. et al. Intron-3 retention/splicing controls neuronal expression of CDTS metrics: http://www.hli-opendata.com/noncoding. apolipoprotein E in the CNS. J. Neurosci. 28, 1452–1459 (2008). phastCons20 metrics: http://hgdownload.cse.ucsc.edu/goldenPath/hg38/ 24. Dieter, L. S. & Estus, S. Isoform of APOE with retained intron 3; quantitation phastCons20way/. and identification of an associated single nucleotide polymorphism. Mol. ROSMAP studies BAM files: https://www.synapse.org/#!Synapse:syn4164376. Neurodegener. 5, 34–34 (2010). OMIM API: http://api.omim.org. 25. GTEx Consortiums. Human genomics. The Genotype-Tissue Expression STOPGAP database: https://github.com/StatGenPRD/STOPGAP/blob/master/ (GTEx) pilot analysis: multitissue gene regulation in humans. Science 348, STOPGAP_data/stopgap.bestld.RData. 648–660 (2015). GTEx portal: https://www.gtexportal.org/home/datasets. 26. Zhang, D. et al. Incomplete annotation has a disproportionate impact on our Ensembl v92: https://www.ensembl.org/index.html. understanding of Mendelian and complex neurogenetic disorders. Sci Adv. 6, GENCODE: https://www.gencodegenes.org/pages/data_access.html. eaay8299 (2020). ENCODE list of problematic regions: https://github.com/Boyle-Lab/. 27. Collado-Torres, L. et al. Flexible expressed region analysis for RNA-seq with Chimpanzee and human bulk RNA-sequencing data: NCBI Gene Expression derfinder. Nucleic Acids Res. 45, e9 (2017). Omnibus; https://www.ncbi.nlm.nih.gov/geo/, accession number GSE127898). 28. Hindorff, L. A. et al. Potential etiologic and functional implications of Source Data are provided with this paper. genome-wide association loci for human diseases and traits. Proc. Natl Acad. Sci. USA 106, 9362–9367 (2009). 29. Savage, J. E. et al. Genome-wide association meta-analysis in 269,867 Received: 28 April 2020; Accepted: 3 March 2021; individuals identifies new genetic and functional links to intelligence. Nat. Genet. 50, 912–919 (2018). 30. Sarropoulos, I., Marin, R., Cardoso-Moreira, M. & Kaessmann, H. Developmental dynamics of lncRNAs across mammalian organs and species. Nature 571, 510–514 (2019). 31. Salta, E. & De Strooper, B. Noncoding RNAs in neurodegeneration. Nat. Rev. References Neurosci. 18, 627–640 (2017). 1. Walker, L. C. & Jucker, M. The exceptional vulnerability of humans to 32. Cardoso-Moreira, M. et al. Gene expression across mammalian organ Alzheimer’s disease. Trends Mol. Med. 23, 534–545 (2017). development. Nature 571, 505–509 (2019). 10 NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 ARTICLE 33. Rathelot, J. A. & Strick, P. L. Subdivisions of primary motor cortex based on (MRC) through the award of a Tenure Track Clinician Scientist Fellowship (MR/ cortico-motoneuronal cells. Proc. Natl Acad. Sci. USA 106, 918–923 (2009). N008324/1). J.H. was supported by the UK Dementia Research Institute which receives 34. de Noordhout, A. M. et al. Corticomotoneuronal synaptic connections in its funding from DRI Limited, funded by the UK Medical Research Council, Alzheimer’s normal man: an electrophysiological study. Brain 122, 1327–1340 (1999). Society and Alzheimer’s Research UK. J.H. has also been funded by the Medical Research 35. Al-Chalabi, A. et al. The genetics and neuropathology of amyotrophic lateral Council (award MR/N026004/1), Wellcome Trust (award 202903/Z/16/Z), Dolby Family sclerosis. Acta Neuropathol. 124, 339–352 (2012). Fund and National Institute for Health Research University College London Hospitals 36. Yu, J. T., Tan, L. & Hardy, J. Apolipoprotein E in Alzheimer’s disease: an Biomedical Research Centre. J.B. is supported through the Science and Technology update. Annu. Rev. Neurosci. 37, 79–100 (2014). Agency, Séneca Foundation, CARM, Spain (research project 00007/COVI/20). The 37. Adusumalli, S., Ngian, Z.-K., Lin, W.-Q., Benoukraf, T. & Ong, C.-T. authors are grateful to Dr Julia di Iulio for providing genomic territory annotation files Increased intron retention is a post-transcriptional signature associated with used to generate Fig. 2a. progressive aging and Alzheimer’s disease. Aging Cell 18, e12928 (2019). 38. Finucane, H. K. et al. Partitioning heritability by functional annotation using genome-wide association summary statistics. Nat. Genet. 47, 1228–1235 Author contributions Z.C. and D.Z. generated the annotation and conducted further analyses. Z.C. and R.H.R. (2015). performed LDSC analysis. Z.C. and E.K.G. generated cDNA and completed Sanger 39. Siepel, A. et al. Evolutionarily conserved elements in vertebrate, insect, worm, sequencing. Z.C. conducted the RNA-sequencing data analyses of APOE. D.Z. and S.G.R. and yeast genomes. Genome Res. 15, 1034–1050 (2005). developed the vizER platform for visualisation of CNC scores. K.D’S., A.F.-B. and J.V. 40. Amemiya, H. M., Kundaje, A. & Boyle, A. P. The ENCODE blacklist: helped with further analyses of RNA-sequencing data. IPDGC contributed PD GWAS identification of problematic regions of the genome. Sci. Rep. 9, 9354 (2019). data summary statistics. J.B., S.G.A.T., H.H. and J.H. provided further guidance on 41. Harrow, J. et al. GENCODE: the reference human genome annotation for the technical aspects of the study. M.R. conceived and supervised the study. Z.C. and M.R. ENCODE project. Genome Res. 22, 1760–1774 (2012). drafted the first version of the manuscript and all authors contributed to the writing and 42. Aken, B. L. et al. Ensembl 2017. Nucleic Acids Res. 45, D635–D642 (2017). reviewing of the manuscript. 43. Bulik-Sullivan, B. K. et al. LD score regression distinguishes confounding from polygenicity in genome-wide association studies. Nat. Genet. 47, 291–295 (2015). Competing interests 44. Jansen, I. E. et al. Genome-wide meta-analysis identifies new loci and The authors declare no competing interests. functional pathways influencing Alzheimer’s disease risk. Nat. Genet. 51, 404–413 (2019). 45. Nalls, M. A. et al. Identification of novel risk loci, causal insights, and heritable Additional information risk for Parkinson’s disease: a meta-analysis of genome-wide association Supplementary information The online version contains supplementary material studies. Lancet Neurol. 18, 1091–1102 (2019). available at https://doi.org/10.1038/s41467-021-22262-5. 46. Pardinas, A. F. et al. Common schizophrenia alleles are enriched in mutation- intolerant genes and in regions under strong background selection. Nat. Correspondence and requests for materials should be addressed to M.R. Genet. 50, 381–389 (2018). Peer review information Nature Communications thanks Terry Goldberg and the other, 47. Wray, N. R. et al. Genome-wide association analyses identify 44 risk variants anonymous, reviewer(s) for their contribution to the peer review of this work. Peer and refine the genetic architecture of major depression. Nat. Genet. 50, reviewer reports are available. 668–681 (2018). 48. Reimand, J. et al. g:Profiler—a web server for functional interpretation of gene Reprints and permission information is available at http://www.nature.com/reprints lists (2016 update). Nucleic Acids Res. 44, W83–W89 (2016). 49. Supek, F., Bosnjak, M., Skunca, N. & Smuc, T. REVIGO summarizes and visualizes long lists of gene ontology terms. PLoS ONE 6, e21800 (2011). Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in 50. Khrameeva, E. et al. Single-cell-resolution transcriptome map of human, published maps and institutional affiliations. chimpanzee, bonobo, and macaque brains. Genome Res. 30, 776–789 (2020). 51. Dobin, A. & Gingeras, T. R. Mapping RNA-seq Reads with STAR. Curr. Protoc. Bioinform. 51, 11.14.11–11.14.19 (2015). Open Access This article is licensed under a Creative Commons 52. Robinson, J. T. et al. Integrative genomics viewer. Nat. Biotechnol. 29, 24–26 Attribution 4.0 International License, which permits use, sharing, (2011). adaptation, distribution and reproduction in any medium or format, as long as you give 53. Bennett, D. A., Schneider, J. A., Arvanitakis, Z. & Wilson, R. S. Overview and appropriate credit to the original author(s) and the source, provide a link to the Creative findings from the religious orders study. Curr. Alzheimer Res. 9, 628–645 Commons license, and indicate if changes were made. The images or other third party (2012). material in this article are included in the article’s Creative Commons license, unless 54. Collado-Torres, L. et al. Reproducible RNA-seq analysis using recount2. Nat. indicated otherwise in a credit line to the material. If material is not included in the Biotechnol. 35, 319–321 (2017). article’s Creative Commons license and your intended use is not permitted by statutory 55. Ng, B. et al. An xQTL map integrates the genetic architecture of the human regulation or exceeds the permitted use, you will need to obtain permission directly from brain’s transcriptome and epigenome. Nat. Neurosci. 20, 1418–1426 (2017). the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Acknowledgements The authors are grateful to the participants in the Religious Order Study, the Memory © The Author(s) 2021 and Aging Project. Z.C. and R.H.R. were supported by grants from the Leonard Wolfson Foundation. M.R. was supported by the United Kingdom Medical Research Council International Parkinson’s Disease Genomics Consortium (IPDGC) Alastair J. Noyce12,13, Rauan Kaiyrzhanov13, Ben Middlehurst14, Demis A. Kia13, Manuela Tan15, , Henry Houlden 8, Huw R. Morris15, Helene Plun-Favreau13, Peter Holmans16, , John Hardy1,4,5,6,7, Daniah Trabzuni13,17, Jose Bras3,13, John Quinn14, Kin Y. Mok13, Kerri J. Kinghorn18, Kimberley Billingsley14, Nicholas W. Wood13, Patrick Lewis19, Sebastian Schreglmann13, Rita Guerreiro13,20, Ruth Lovering21, Lea R’Bibo13, Claudia Manzoni19, Mie Rizig13, & Mina Ryten 1,2,3 ✉, Sebastian Guelfi1, Valentina Escott-Price22, Viorica Chelban13, Thomas Foltynie15, Nigel Williams22, Alexis Brice23, Fabrice Danjou23, Suzanne Lesage23, NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 Jean-Christophe Corvol24, Maria Martinez25, Claudia Schulte26, Kathrin Brockmann26, Javier Simón-Sánchez26, Peter Heutink26, Patrizia Rizzu26, Manu Sharma27, Thomas Gasser26, Aude Nicolas28, Mark R. Cookson28, Sara Bandres-Ciga28, Cornelis Blauwendraat28,29, David W. Craig30, Faraz Faghri28,31, J. Raphael Gibbs28, Dena G. Hernandez28, Kendall Van Keuren-Jensen32, Joshua M. Shulman33,34, Hampton L. Leonard28, Mike A. Nalls28,35, Laurie Robak33,34, Steven Lubbe36, Steven Finkbeiner37,38,39, Niccolo E. Mencacci40, Codrin Lungu41, Andrew B. Singleton28, Sonja W. Scholz42, Xylena Reed28, Roy N. Alcalay43,44, Ziv Gan-Or45, Guy A. Rouleau45, Lynne Krohn45, Jacobus J. van Hilten46, Johan Marinus46, Astrid D. Adarmes-Gómez47, Miquel Aguilar48, Ignacio Alvarez48, Victoria Alvarez49, Francisco Javier Barrero50, Jesús Alberto Bergareche Yarza51, Inmaculada Bernal-Bernal52, Marta Blazquez49, Marta Bonilla-Toribio52, Juan A. Botía8, María Teresa Boungiorno48, Dolores Buiza-Rueda47, Ana Cámara52, Fátima Carrillo47, Mario Carrión-Claro47, Debora Cerdan53, Jordi Clarimón54,55, Yaroslau Compta52, Monica Diez-Fairen48, Oriol Dols-Icardo54,55, Jacinto Duarte53, Raquel Duran56, Francisco Escamilla-Sevilla57, Mario Ezquerra52, Cici Feliz58, Manel Fernández52, Rubén Fernández-Santiago52, Ciara Garcia49, Pedro García-Ruiz58, Pilar Gómez-Garre47, Maria Jose Gomez Heredia59, Isabel Gonzalez-Aramburu60, Ana Gorostidi Pagola51, Janet Hoenicka61, Jon Infante55,62, Silvia Jesús47, Adriano Jimenez-Escrig62, Jaime Kulisevsky55,63, Miguel A. Labrador-Espinosa47, Jose Luis Lopez-Sendon62, Adolfo López de Munain Arregui51, Daniel Macias47, Irene Martínez Torres64, Juan Marín54,55, Maria Jose Marti52, Juan Carlos Martínez-Castrillo62, Carlota Méndez-del-Barrio47, Manuel Menéndez González49, Marina Mata65, Adolfo Mínguez57, Pablo Mir47, Elisabet Mondragon Rezola51, Esteban Muñoz47, Javier Pagonabarraga63, Pau Pastor48, Francisco Perez Errazquin59, Teresa Periñán-Tocino47, Javier Ruiz-Martínez66, Clara Ruz56, Antonio Sanchez Rodriguez60, María Sierra60, Esther Suarez-Sanmartin49, Cesar Tabernero53, Juan Pablo Tartari67, Cristina Tejera-Parrado47, Eduard Tolosa52, Francesc Valldeoriola52, Laura Vargas-González47, Lydia Vela68, Francisco Vives56, Alexander Zimprich69, Lasse Pihlstrom70, Mathias Toft70, Sulev Koks71,72, Pille Taba73 & Sharon Hassin-Baer74,75 12 Preventive Neurology Unit, Wolfson Institute of Preventive Medicine, Barts and the London School of Medicine and Dentistry, Queen Mary University of London, London, UK. 13Department of Molecular Neuroscience, Queen Square Institute of Neurology, University College London (UCL), London, UK. 14Institute of Translational Medicine, University of Liverpool, Liverpool, UK. 15Department of Clinical and Movement Neuroscience, Queen Square Institute of Neurology, University College London (UCL), London, UK. 16Biostatistics & Bioinformatics Unit, Institute of Psychological Medicine and Clinical Neuroscience, MRC Centre for Neuropsychiatric Genetics & Genomics, Cardiff, UK. 17Department of Genetics, King Faisal Specialist Hospital and Research Centre, Riyadh, Saudi Arabia. 18Institute of Healthy Ageing, University College London (UCL), London, UK. 19University of Reading, Reading, UK. 20UK Dementia Research Institute, University College London (UCL), London, UK. 21 Institute of Cardiovascular Science, University College London (UCL), London, UK. 22MRC Centre for Neuropsychiatric Genetics and Genomics, Cardiff University School of Medicine, Cardiff, UK. 23Institut 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. 24Institut 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, Centre d’Investigation Clinique Pitié Neurosciences CIC-1422, AP-HP, Pitié-Salpêtrière Hospital, Paris, France. 25INSERM UMR 1220 and Paul Sabatier University, Toulouse, France. 26 Department for Neurodegenerative Diseases, Hertie Institute for Clinical Brain Research, German Center for Neurodegenerative Diseases, DZNE, German Center for Neurodegenerative Diseases, University of Tübingen, Tübingen, Germany. 27Centre for Genetic Epidemiology, Institute for Clinical Epidemiology and Applied Biometry, University of Tubingen, Tübingen, Germany. 28Laboratory of Neurogenetics, National Institute on Aging, Bethesda, MD, USA. 29National Institute of Neurological Disorders and Stroke, Bethesda, MD, USA. 30Department of Translational Genomics, Keck School of Medicine, University of Southern California, Los Angeles, CA, USA. 31Department of Computer Science, University of Illinois at Urbana-Champaign, Urbana, IL, USA. 32Neurogenomics Division, TGen, Phoenix, AZ, USA. 33Departments of Neurology, Neuroscience, and Molecular & Human Genetics, Baylor College of Medicine, Houston, TX, USA. 34Jan and Dan Duncan Neurological Research Institute, Texas Children’s Hospital, Houston, TX, USA. 35Data Tecnica International, Glen Echo, MD, USA. 36Ken and Ruth Davee Department of Neurology, Northwestern University Feinberg School of Medicine, Chicago, IL, USA. 37Department of Neurology and Physiology, University of California, San Francisco, CA, USA. 38Gladstone Institute of Neurological Disease, San Francisco, CA, USA. 39Taube/Koret Center for Neurodegenerative Disease Research, San Francisco, CA, USA. 40Northwestern University Feinberg School of Medicine, Chicago, IL, USA. 41National Institutes of Health Division of Clinical Research, NINDS, National Institutes of Health, Bethesda, MD, USA. 42Neurodegenerative Diseases Research Unit, National Institute of Neurological Disorders and Stroke, Bethesda, MD, USA. 43Department of Neurology, College of Physicians and Surgeons, Columbia University Medical Center, New York, NY, USA. 44Taub Institute for Research on Alzheimer’s Disease and the Aging Brain, College of Physicians and Surgeons, Columbia University Medical Center, New York, NY, USA. 45Montreal Neurological Institute and Hospital, Department of Neurology & Neurosurgery, Department of Human Genetics, McGill University, Montréal, QC, Canada. 46Department of Neurology, Leiden University Medical Center, Leiden, The Netherlands. 47Instituto de Biomedicina de Sevilla Hospital Universitario Virgen del Rocío/CSIC/Universidad de 12 NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-021-22262-5 ARTICLE Sevilla, Seville, Spain. 48Fundació Docència I Recerca Mútua de Terrassa and Movement Disorders Unit, Department of Neurology, University Hospital Mutua de Terrassa, Terrassa, Barcelona, Spain. 49Hospital Universitario Central de Asturias, Oviedo, Spain. 50Hospital Universitario San Cecilio de Granada, Universidad de Granada, Granada, Spain. 51Instituto de Investigación Sanitaria Biodonostia, San Sebastián, Spain. 52Hospital Clinic de Barcelona, Barcelona, Spain. 53Hospital General de Segovia, Segovia, Spain. 54Memory Unit, Department of Neurology, IIB Sant Pau, Hospital de la Santa Creu i Sant Pau, Universitat Autònoma de Barcelona, Barcelona, Spain. 55Centro de Investigación Biomédica en Red en Enfermedades Neurodegenerativas, Madrid, Spain. 56Centro de Investigacion Biomedica, Universidad de Granada, Granada, Spain. 57Hospital Universitario Virgen de las Nieves, Instituto de Investigación Biosanitaria de Granada, Granada, Spain. 58Hospital Universitario Marqués de Valdecilla-IDIVALS, Instituto de Investigación Sanitaria Fundación Jiménez Díaz, Madrid, Spain. 59Hospital Universitario Virgen de la Victoria, Malaga, Spain. 60Hospital Universitario Marqués de Valdecilla-IDIVAL, Santander, Spain. 61Institut de Recerca Sant Joan de Déu, Terrassa, Barcelona, Spain. 62Hospital Universitario Ramón y Cajal, Madrid, Spain. 63Movement Disorders Unit, Department of Neurology, IIB Sant Pau, Hospital de la Santa Creu i Sant Pau, Universitat Autònoma de Barcelona, Terrassa, Barcelona, Spain. 64Department of Neurology, Instituto de Investigación Sanitaria La Fe, Hospital Universitario y Politécnico La Fe, Valencia, Spain. 65Department of Neurology, Hospital Universitario Infanta Sofía, Madrid, Spain. 66Hospital Universitario Donostia, Instituto de Investigación Sanitaria Biodonostia, San Sebastián, Spain. 67Fundació Docència i Recerca Mútua de Terrassa and Movement Disorders Unit, Department of Neurology, University Hospital Mutua de Terrassa, Terrassa, Barcelona, Spain. 68Department of Neurology, Hospital Universitario Fundación Alcorcón, Madrid, Spain. 69Department of Neurology, Medical University of Vienna, Vienna, Austria. 70Department of Neurology, Oslo University Hospital, Oslo, Norway. 71Department of Pathophysiology, University of Tartu, Tartu, Estonia. 72Department of Reproductive Biology, Estonian University of Life Sciences, Tartu, Estonia. 73Department of Neurology and Neurosurgery, University of Tartu, Tartu, Estonia. 74Movement Disorders Institute, Department of Neurology and Sagol Neuroscience Center, Chaim Sheba Medical Center, Ramat Gan, Israel. 75Sackler Faculty of Medicine, Tel Aviv University, Tel Aviv, Israel. NATURE COMMUNICATIONS | (2021)12:2076 | https://doi.org/10.1038/s41467-021-22262-5 | www.nature.com/naturecommunications 13