Whole exome sequencing study identifies novel rare and common Alzheimer’s-Associated variants involved in immune response and transcriptional regulation

Molecular Psychiatry https://doi.org/10.1038/s41380-018-0112-7 ARTICLE Whole exome sequencing study identifies novel rare and common Alzheimer’s-Associated variants involved in immune response and transcriptional regulation Joshua C. Bis1et al ●Alzheimer’s Disease Sequencing Project Received: 21 December 2017 / Revised: 1 May 2018 / Accepted: 14 May 2018 © The Author(s) 2018. This article is published with open access Abstract The Alzheimer’s Disease Sequencing Project (ADSP) undertook whole exome sequencing in 5,740 late-onset Alzheimer disease (AD) cases and 5,096 cognitively normal controls primarily of European ancestry (EA), among whom 218 cases and 177 controls were Caribbean Hispanic (CH). An age-, sexand APOE based risk score and family history were used to select cases most likely to harbor novel AD risk variants and controls least likely to develop AD by age 85 years. We tested ~1.5 million single nucleotide variants (SNVs) and 50,000 insertion-deletion polymorphisms (indels) for association to AD, using multiple models considering individual variants as well as gene-based tests aggregating rare, predicted functional, and loss of function variants. Sixteen single variants and 19 genes that met criteria for significant or suggestive associations after multiple-testing correction were evaluated for replication in four independent samples; three with whole exome sequencing (2,778 cases, 7,262 controls) and one with genome-wide genotyping imputed to the Haplotype Reference Consortium panel (9,343 cases, 11,527 controls). The top findings in the discovery sample were also followed-up in the ADSP whole-genome sequenced family-based dataset (197 members of 42 EA families and 501 members of 157 CH families). We identified novel and predicted functional genetic variants in genes previously associated with AD. We also detected associations in three novel genes: IGHG3 (p =9.8 × 10−7), an immunoglobulin gene whose antibodies interact with β-amyloid, a long noncoding RNA AC099552.4 (p =1.2 × 10−7), and a zinc-finger protein ZNF655 (gene-based p =5.0 × 10−6). The latter two suggest an important role for transcriptional regulation in AD pathogenesis. Introduction Genomic studies have revealed that late-onset Alzheimer disease (LOAD) is highly polygenic, with as many as 30 susceptibility loci identified through large-scale metaanalysis of genome-wide association studies (GWAS), targeted exome genotyping array, and several early whole exome sequencing (WES) studies [1–12]. Although AD susceptibility is highly heritable (h2=0.58–0.79) [13], much of its genetic architecture is still unknown and few rare variants have been detected thus far [3,6,7,14–19]. Discovery of rare variants in genomic studies, even those with large sample sizes and examining highly heritable diseases, remains challenging due to statistical power limitations in detecting all but the most strongly associated variants (odds ratio (OR) > 1.5) [20–23]. The protein coding regions of the genome, or exome, are the best characterized and most conserved portions of the genome and the source of most variants identified to date that are responsible for Mendelian diseases [24]; thus, the exome is a more Alzheimer’s Disease Sequencing Project members are listed below the Acknowledgement These authors contributed equally: Joshua C. Bis, Xueqiu Jian, Brian W. Kunkle, Yuning Chen These authors equally supervised the study: Adam C. Naj, Myriam Fornage, Lindsay A. Farrer *Lindsay A. Farrer [email protected] Extended author information available on the last page of the article Electronic supplementary material The online version of this article (https://doi.org/10.1038/s41380-018-0112-7) contains supplementary material, which is available to authorized users. Variant summary data can be found at the NIA Genetics of Alzheimer 's Disease Data Storage site (https://www.niagads.org) under accession number NG00065. 1234567890();,: 1234567890();,: attractive and less expensive target for identifying rare variants of large effect on disease than the non-protein coding portion of the genome. The Alzheimer’s Disease Sequencing Project (ADSP) was developed jointly by the National Institute on Aging (NIA) and National Human Genome Research Institute (NHGRI) in response to the National Alzheimer’s Project Act milestones (https://aspe.hhs.gov/national-alzheimers-project-act)tofight Alzheimer’s disease (AD) as an effort to analyze the genomes of well-characterized individuals with and without AD. To detect rare variants and genes associated with LOAD, we performed single-variant and gene-based analyses, including annotated loss-of-function analyses, on the ADSP Discovery Phase Case-Control WES dataset, and attempted to replicate associations in three independent WES datasets, a GWAS dataset containing single nucleotide variants (SNVs) that were imputed using the Haplotype Reference Consortium (HRC) [25] reference panel, and the ADSP family-based whole genome sequence dataset. Methods Sample selection and data preparation Study participants were either European-American (EA) or Caribbean Hispanic (CH) ancestry and were sampled in two ways. To maximize contrast between cases and controls, and power to discover novel associations, the majority of participants were chosen using a risk score that included dosages of the APOE ε2/ε3/ε4 alleles, sex and either onset age (for cases) or age at last exam for controls (or pathology-based adjusted age at death for neuropathology control) [26]. All cases were at least 60 years old and met NINCDS-ADRDA criteria for possible, probable or definite AD based on clinical assessment, or had presence of AD (moderate or high likelihood) upon neuropathology examination. To maximize our ability to discover novel genetic associations, we chose cases whose AD risk score indicated that their disease was not well explained by age, sex, or dosages of the APOE ε2/ε3/ε4 alleles. Conversely, cognitively healthy controls were selected with the goal of identifying alleles associated with the increased risk of or protection from lateonset AD. At the time of last exam, all potential controls were at least 60 years old and were either judged to be cognitively normal or did not meet pathological criteria [27,28] for AD following brain autopsy. Controls were selected for this study on the basis of the risk score indicating that they were the least likely to develop AD by age 85 years. Applying the risk score resulted in a sample that contained 2,220 AD cases (40%) and 752 controls (14%) who were ε4 heterozygotes and 161 AD cases (3%) and 17 controls (<1%) who were ε4 homozygotes. In addition, we sampled a set of “enriched”cases from families having at least three affected members for whom the diagnosis of AD was verified by direct examination or review of cognitive testing data and medical records. Cases from early-onset AD families or families with a known PSEN1,PSEN2,orAPP mutation were excluded. Within each family, we selected only one AD case, typically the member with the lowest a priori AD risk (based on the risk score defined above), provided this person had sufficient genomic DNA. In addition, because 172 of the “enriched” cases described above were of CH ancestry, we also selected a set of 171 ageand sex-matched cognitively normal CH participants to serve as controls. Participant characteristics are shown in Table 1A. Procedures Genotype calling and data processing WES data were generated at the Broad Institute, the Baylor College of Medicine’s Human Genome Sequencing Center, and Washington University’s McDonnell Genome Institute. An effort was made to assign all samples from a study of origin to the same center and there was a relatively balanced number of cases and controls at each center. Genotypes for bi-allelic SNVs and insertion-deletion polymorphisms (indels) were called using ATLAS2 using version Ch37/ hg19 of the reference genome. A coordinated effort was implemented for centralized variant calling and quality control (QC) efforts in order to create one batch of data for analysis. Although there were differences in allele frequencies across sequencing centers for some variants, it was difficult to determine whether these represented technical artifacts of different capture kits, variability in genetic background among cohorts assembled for this study, or chance differences that will often occur for infrequent or rare variants. QC steps and methods for evaluating cryptic relatedness, population substructure, differential missingness, and variant annotation are described in the Supplementary Materials. Single-variant and gene-based association analyses Statistical models & rationale for covariate adjustments All models included adjustment for sequencing center and population substructure. Before conducting the primary analyses, we evaluated up to 20 PCs for association with AD status. Only ancestry-specific PCs that significantly associated with AD status (P<0.005) in at least one of the three adjustment models were included as covariates (EA subgroup: PC1, PC5, PC8, PC9, PC10, PC11, PC18; J. C. Bis et al. Hispanic analyses included PC1 and PC2). Because most participants for the discovery study were sampled to maximize differences in cases and controls based on age, sex, and APOE genotypes, we included only PCs and sequencing center in our base adjustment model (Model 0). We evaluated two other models that included several covariates in addition to those in the base model: Model 1 adjusted for sex and age at diagnosis or last follow-up; and Model 2 adjusted for APOE ε4&ε2 dosages in addition to those included in Model 1. All analyses were performed separately by ancestry (EA and CH) using seqMeta (version 1.6) [29]; the primary analysis is an inverse variance-weighted meta-analysis of these two groups. Single variant tests were limited to variants with at least 10 copies of the minor allele across the total QCed sample (MAF~0.0005). Gene-based association testing Gene-based tests examine the aggregate effect of risk and protective variants within a region defined by gene annotations. We performed gene-based tests using SKAT-O, which optimally combines SKAT and burden tests [30]. For these analyses, the SKAT portion of the test included variants with a MAF≤0.05; the burden component aggregated variants with MAF≤0.01. The SKAT test used ‘Wu weights’,defined by a beta density function with prespecified parameters a1 =1 and a2 =25, evaluated at the sample minor allele frequency. The SKAT-O statistic, a linear combination between a SKAT statistic (Qskat) and a burden statistic (Qburden) equal to (1-ρ)Q skat +ρQburden, was optimized across 11 values of ρ(0.1 increments), and calculation of the significance took into consideration the multiple values of ρevaluated. In order to improve power by removing variants predicted to have a low functional impact on the translated protein, we filtered variants in each gene on the basis of annotated function as described in the Supplementary Materials. We performed SKAT-O testing for genes with at least two qualifying variants contributing to the test. The minimum number of aggregated alleles (i.e., cumulative minor allele count or cMAC) for a gene-based test was set at 10. Statistical significance thresholds for discovery stage analyses Within each analysis framework including individual variants and gene-based aggregation of variants evaluated under particular functional annotation criteria, suggestive associations (p < 1/ # tests) were selected for follow-up testing in independent samples and a Bonferroni-corrected threshold was used to define experiment wide statistical significance (p < 0.05 / # tests). We did not correct for the three models and meta-analyses of the combined results of the EA +CH populations because the results were highly correlated across the covariate adjustment strategies (Supplementary Figure S2). Replication sample and analyses Primary replication analyses for the SNVs / genes that we identified to be genome-wide significant or suggestive in any model were conducted in three independent WES Table 1 Participant Characteristics A. Discovery Sample AD Cases (N =5,740) Cognitively Normal Controls (N =5,096) Ancestry Sampling N Age (mean) Sex (%F) APOE E4 (%carrier) N Age (mean) Sex (%F) APOE E4 (%carrier) EA Case-Control 5,015 75.25 55.8% 40.6% 4,919 86.53 59.2% 14.4% EA Enriched 507 83.61 63.3% 66.9% NA NA NA NA Hispanic Case-Control 46 72.59 71.7% 43.5% 6 85.94 66.7% 16.7% Hispanic Enriched 172 75.45 61.6% 39.5% 171 73.46 60.8% 39.2% B. Replication Sample AD Cases (N =12,121) Cognitively Normal Controls (N =18,789) N Age (mean) Sex (% F) APOE ε4 (% carrier) N Age (mean) Sex (% F) APOE ε4 (% carrier) CHARGE WES 612 81 67% 54% 1,836 80 58% 24% ADES-FR WES 1,142 74 64% 49% 1,104 80 58% 22% FinnAD WES 1,024 74 62% -- 4,322 71 51% -- ADGC GWAS 9,343 74 54% 64% 11,527 75 54% 25% Whole exome sequencing study identifies novel rare and common Alzheimer’s-Associated variants. . . datasets including CHARGE (612 cases and 1,836 controls), ADES-FR (1,142 cases and 1,104 controls) [31], and FinnAD (1,024 cases and 4,322 controls), as well as in the Alzheimer’s Disease Genetics Consortium (ADGC) HRCimputed GWAS dataset (Table 1B, Supplementary Materials, Supplementary Table S1). The ADGC dataset included GWAS data on 9,343 cases and 11,527 cognitively normal elders from 32 datasets for whom genotypes were imputed using the Haplotype Reference Consortium (HRC) r1.1 reference panel (Supplementary Table S2) [25,32]. CHARGE and ADGC participants selected for ADSP discovery analyses were not included in the replication study. Because we included all available cases and controls in the replication datasets instead of applying the participant selection criteria used for the discovery sample to maximize difference in cases and controls, model 0 is not appropriate in replication studies. Hence, single variant tests and genebased SKAT-O tests were performed using seqMeta for models 1 and 2 only. Meta-analysis of summarized results from the four samples was performed using seqMeta. We also performed a meta-analysis of results across the ADSP discovery and four replication cohorts for findings that were at least “suggestive”(p < 1 / genes or variants) in the discovery phase. In addition to models 1 and 2, we conducted a meta-analysis of results obtained using model 0 in the ADSP discovery data and model 1 in the replication cohorts to verify our findings in ADSP model 0. Because the ADSP discovery dataset includes CH participants and all replication cohorts consist of EA participants only, we performed the meta-analysis with and without CH participants in ADSP. We considered any variants or genes with two-stage meta-analysis p-values < 0.05 / # tests to be significant per the recommendation by Skol et al. that joint analysis is more efficient than replication-based analysis for two-stage genome-wide association studies [33]. We acknowledge, however, that additional replication in independent samples is required. The top findings in the discovery sample were also followed-up in the ADSP whole-genome sequenced (WGS) family-based dataset [34,35]. This dataset includes 197 individuals sequenced in 42 EA families and 351 individuals in 67 CH families. Additional follow-up was also performed using WGS information from 150 members of another 48 CH families. Individual variants were assessed by examining their co-segregation with AD status within families. Gene-based association tests were performed using the FSKAT software [36]. Analysis of variants at previously established AD loci To identify a set of variants related to AD risk in loci previously associated with AD, we compiled a list of genes containing variants with significant or suggestive associations (p < 1 × 10−3) in either the published IGAP or UKBB AD GWAS meta-analyses [9,37]. Because many signal variants from GWAS are in intergenic regions, we used a combination of BEDOPS [38] and BEDTools [39] operations to enlarge the genomic coordinates of these associated variants by 50 Kb on each side, merging adjacent regions that were overlapping and/or book-ended. Of the resulting genomic regions, segments greater than 100 kb were retained, shortened by 50 kb on each side, merged if separated by 200 kb or less, and utilized to find overlapping protein-coding genes, with gene boundaries as defined in version 19 of the GENCODE gene set [40] and a 50 kb buffer on each side. These parameters and sequence of operations were chosen because they resulted in an algorithm that satisfactorily captured the genomic interval of the association landscape at each locus, as confirmed by visual inspection of LocusZoom regional plots [41]. We queried variant and gene level association statistics for the resulting list of 299 putatively associated AD genes. Results Description of study samples after QC and filtering After exclusions, 10,836 participants were available for analysis (5,740 cases; 5,096 controls). This included 218 CH cases and 177 CH controls. The study included more women than men, and, due largely to the selection criteria, cases were younger on average than controls and were more likely to carry one copy of the APOE ε4 allele. In total, the data included 1,524,414 bi-allelic SNVs or short indels. Most variants were rare, with 1,493,926 (98%) of variants having minor allele frequency of less than 5% and 160,898 (11%) having a minor allele count (MAC) of at least 10 copies. Single-variant SNV and short indels association analysis We performed single variant analyses for the 160,898 variants with a combined minor allele count of at least 10 copies across all participants (Supplementary Figure S3). Genomic inflation was moderate (λ< 1.1 in all models) (Supplementary Figures S4–S7). Single variant association testing identified three variants at an exome-wide significance level (p < 3.1 × 10−7) and 14 variants at the suggestive threshold (p < 6.1 × 10−6) outside of the APOE region (Table 2, Fig. 1, Supplementary Table S3). The significant associations included the rare missense R47H variant in TREM2 (rs75932628, p =4.8 × 10−12), a common variant in PILRA (rs2405442, p =1.7 × 10−7), and a novel rare variant in the long non-coding RNA AC099552.4 J. C. Bis et al. Table 2 Associations with Individual Variants outside the APOE region ADSP Discovery Meta All Replication (ADGC+CHARGE +ADES-FR+FinnAD) Discovery +All Replication Name gene MAC (EA/CH) best P Model (group) MAC P Model 1 P Model 2 P Model 0 P Model 1 P Model 2 6:41129252:C:T (R47H) TREM2 120/0 4.8E-12 0 (EA) 224 1.6E-06 2.7E-06 3.2E-16 2.8E-10 1.6E-10 7:154988675:G:A AC099552.4 10/0 1.2E-07 2 (EA) 0 NA NA 1.3E-02 2.0E-07 1.2E-07 7:99971313:T:C (rs2405442) PILRA 6,219/219 1.7E-07 0 (EA) 22,798 5.3E-05 2.3E-05 9.5E-10 1.1E-06 5.0E-07 20:62729814:C:T (rs148484121) OPRL1 61/4 5.8E-07 1 (all) 111 3.4E-01 5.6E-01 3.7E-03 1.4E-04 4.5E-04 11:59940599:T:A (rs7232) MS4A6A 7,540/258 7.7E-07 0 (all) 20,963 1.4E-11 3.1E-09 5.6E-17 3.8E-14 2.6E-11 17:44828931:G:A (rs199533) NSF 4,238/135 1.3E-06 0 (all) 11,120 2.5E-01 1.4E-02 2.1E-04 1.6E-02 1.9E-04 14:106235767:C:T (rs77307099) IGHG3 6,200/176 1.9E-06 0 (all) 721 4.0E-01 3.5E-01 1.4E-06 1.3E-04 7.9E-05 14:106235766:G:A (rs78376194) IGHG3 6,202/176 1.9E-06 0 (all) 719 4.2E-01 3.6E-01 1.5E-06 1.4E-04 8.5E-05 6:15638035:C:T (rs77460377) DTNBP1 16/3 1.9E-06 2 (all) 35 8.5E-01 8.7E-01 8.7E-02 5.2E-03 3.0E-03 6:33041297:G:A (rs112178281) HLA-DPA1 10/0 2.9E-06 1 (EA) 6 7.5E-01 9.2E-01 1.4E-01 2.1E-05 2.0E-05 11:59945745:T:C (rs12453) MS4A6A 8,265/258 3.2E-06 0 (EA) 23,420 4.7E-11 3.0E-08 1.4E-15 6.0E-13 1.2E-09 3:195506101:T:A MUC4 38/6 3.8E-06 1 (all) 0 NA NA 3.0E-04 3.8E-06 8.1E-06 10:88446985:T:C (rs76615432) LDB3 760/62 5.0E-06 1 (CH) 2,303 5.3E-01 5.9E-01 7.0E-01 6.3E-01 6.4E-01 19:1047507:AGGAGCAG:A ABCA7 67/0 4.3E-06 0 (EA) 11 8.8E-02 9.7E-02 2.4E-04 1.6E-02 1.7E-02 14:106236128:T:A (rs12890612) IGHG3 6,395/369 4.5E-06 0 (all) 1,473 8.5E-02 7.5E-02 9.8E-07 8.0E-05 6.4E-05 7:99799845:T:A (rs104395) STAG3 5,248/248 5.5E-06 0 (EA) 15,948 3.0E-03 1.2E-03 8.8E-07 1.2E-04 4.0E-05 Table shows variants with P< 6.1 × 10−6in EA, CH,or combined strata in the discovery sample. Exome-wide significant results (P< 3.1 × 10−7) and suggestive results which improved in meta analysis of discovery +replication data are highlighted in bold. Results without variation data in the replication datasets are indicated in italics Whole exome sequencing study identifies novel rare and common Alzheimer’s-Associated variants. . . (7:154,988,675, p =1.2 × 10−7). These results were attenuated when including age, sex, and APOE ε2/ε3/ε4 allele as covariates. Gene-based association analysis combining SNVs and indels We aggregated 918,053 variants with a combined MAF < 0.05 and annotated as high or moderate impact into genebased tests using SKAT-O. This corresponds to 17,613 genes with more than one variant and a cumulative minor allele count (cMAC) of at least 10 copies. Applying more stringent filtering, we limited to variants annotated as high impact; aggregating 42,502 rare or uncommon (MAF < 0.05) variants into 4,634 genes (again, limiting to genes with >1 variant and a cMAC ≥10). For the purposes of identifying novel associations, we considered all genes or variants within 250kb of APOE as part of the APOE locus. Three known genes (ABCA7,TREM2 and CBLC in the APOE region) and two novel genes (OPRL1 and GAS2L2) achieved exome-wide statistical significance for their respective tests in the discovery analyses (Table 3, Fig. 2). Four additional genes (ZNF655,RHBDD1,SIRPB1, and RPS16) reached suggestive significance across the nine models (Fig. 2, Supplementary Table S4). Analyses filtered to include only variants with CADD scores ≥15 or ≥20 produced most of the same top-ranked results as the VEP gene-based results (Supplementary Table S5), noting that the overall VEP High/Moderate and CADD≥15 results, as well as the VEP High and CADD≥20 results, are only moderately correlated (Spearman rank correlation r =0.51) (Supplementary Figure S8). Three novel genes (CACNB3, HHIP-AS1, and RP11-68L1.1) were exome-wide statistically significant in the CH group in analyses restricted to variants with CADD scores ≥20 (Supplementary Table S5), however these are likely false positives because in each instance the result is accounted for by a single variant that was observed in one person only. Loss-of-function (LOF) association analysis Among 78,529 unfiltered high impact variants, 72,694 were annotated as LoF and 68,121 were further deemed as highconfidence, most of which were frameshift and stop-gained (Fig. 3). As expected, over 90% of these high-confidence LoF variants were singletons (53,120, 78%), doubletons (6,579, 10%), or tripletons (2,222, 3%), and most of these were observed in European Americans only. Association analysis of 2,378 high-confidence LoF variants with MAC ≥ 10 with adjustment for sequencing center and PCs revealed one Bonferroni corrected significant p < 2.1 × 10−5) variant, a previously reported frameshift deletion in ABCA7 (Table 3)[42]. Gene-based analysis of 32,863 highconfidence LoF variants with MAF ≤5% mapping to 3,558 genes with at least two variants and cMAC ≥10 (Supplementary Table S4) also showed that ABCA7 with adjustment for sequencing center and PCs (Model 0) and GAS2L2 with adjustment for sequencing center, PCs, sex, age, and APOE genotype (Model 2) reached experimentwide significance threshold p < 1.4 × 10−5). Replication analysis Of the 16 single variants outside the APOE region tested in the replication samples (Table 2, Supplementary Table S3), the TREM2 R47H mutation and four variants in three other previously known genes (one missense and one synonymous variant in MS4A6A, a synonymous variant in PILRA, Fig. 1 Manhattan plot showing genome-wide association results for individual common variants. The plot shows the p-values from the Discovery meta-analysis against their genomic position for association with AD. Only variants with a combined minor allele count of ≥10 were included; the minimum p-value from the three adjustment models for either the meta-analysis, European Ancestry (EA), or Caribbean Hispanic (CH) is plotted for each variant. Genes containing the variant are indicated above points that surpassed our significance threshold for follow-up. The dotted line indicates the threshold for follow-up, p < 6.1 × 10−6, corresponding to (1 / #variants) tested. The dashed line indicates the threshold for exome-wide significance, p < 3.1 × 10−7, corresponding to (0.05 / #variants tested) J. C. Bis et al. and a missense variant in CR1) were significantly associated with AD in the combined discovery and replication analysis (Table 2). Associations with two variants in a novel gene STAG3 (rs1043915, p =5.5 × 10−6) were also replicated and significantly associated with AD. We were unable to assess replication with the novel AC099552.4 variant because it was not observed or imputed in the replication datasets. One of the IGHG3 variants (rs12890621) showed borderline evidence for association in models 1 and 2 (p = 0.085 and 0.075, respectively), and evidence for association was strengthened to near exome-wide significance (p = 9.8 × 10−7) in the combined discovery and replication sample. In total, 19 genes across the nine models were significantly or suggestively associated with AD and were tested in the replication stage (Supplementary Table S4). Gene-based tests including high or moderate impact variants showed evidence for replication and reached genomewide significance in the combined discovery and replication analysis for three genes:TREM2 and two genes in the APOE region (CBLC and BCAM) (Table 3, Supplementary Table S6). The association with GAS2L2, the potential novel gene identified in a model 2 SKAT-O test with high impact SNVs in the discovery sample, was slightly above the nominal significance level (p =0.051 and p =0.067, respectively, in models 1 and 2) in meta-analysis across four replication cohorts. However, this association was only nominally significant in the meta-analysis combining the discovery and replication cohorts (p =0.029 for model 2). In gene-based tests including only high impact SNVs, the known AD risk gene ABCA7 and the potential novel gene ZNF655 reached the nominal p-value of 0.05 in metaanalysis with replication cohorts as well as in a metaanalysis of discovery and replication samples (Table 3, Supplementary Table S6). These two genes were also nominally significant in SKAT-O tests, limited to high impact variants (ZNF655:p=7.9 × 10−6;ABCA7: p=6.2 × 10−5) and LOF variants (ZNF655:p=5.0 × 10−6). Because PILRA, a previously established AD gene [43] is proximate to STAG3 (159 kb) and ZNF655 (797 kb), we performed conditional analysis in the discovery sample to determine whether these novel association signals are independent. These analyses demonstrated that the association with multiple rare variants ZNF655 in the genebased test is distinct from those with common variants in PILRA (p =1.08 × 10−4) and STAG3 (p =7.75 × 10−5). In a model containing both PILRA and STAG3 variants the association with PILRA remains significant (p =0.011) but the association with STAG3 does not (p =0.21) (Supplementary Table S7). Follow-up in ADSP family-based data We also followed-up the significant and suggestive singlevariant and gene-based results from the discovery stage in the Table 3 Gene-based Association Results ADSP Discovery Meta All Replication Discovery +All Replication Gene Variants b SNVs best P Model SNVs P Model 1 P Model 2 SNVs P Model 0 P Model 1 P Model 2 TREM2 High-Mod 50 1.8E-11 0339.3E-10 5.4E-09 65 2.0E-17 3.8E-11 6.0E-11 CBLC aHigh-Mod 44 1.1E-07 0352.5E-20 6.7E-03 61 1.0E-27 6.1E-22 4.9E-02 OPRL1 High-Mod 42 2.6E-06 1 37 1.3E-01 3.0E-01 64 8.3E-03 5.4E-04 1.7E-03 CBX3 High-Mod 8 6.0E-05 0 10 1.3E-01 2.8E-01 17 4.9E-04 4.6E-02 6.1E-02 BCAM aHigh-Mod 90 5.2E-04 1 88 4.7E-19 3.7E-03 144 3.5E-27 2.8E-20 4.8E-02 GAS2L2 High 7 3.9E-06 2 5 5.1E-02 6.7E-02 10 4.5E-01 3.9E-02 2.9E-02 ZNF655 High 9 2.8E-05 0 6 3.2E-02 3.4E-02 13 7.9E-06 8.4E-04 3.4E-04 RHBDD1 High 2 3.2E-05 2 4 8.8E-01 9.8E-01 5 3.5E-01 4.8E-01 2.7E-01 SIRPB1 High 6 8.0E-05 2 3 9.2E-01 7.9E-01 8 6.4E-01 3.0E-01 2.6E-01 RPS16 High 5 1.6E-04 2 2 7.4E-01 4.2E-01 5 4.4E-02 7.7E-03 6.5E-03 ABCA7 LoF 43 2.1E-06 0 16 1.5E-01 1.1E-01 51 1.2E-04 1.2E-03 3.4E-04 GAS2L2 LoF 7 3.9E-06 2 3 3.9E-02 4.8E-02 8 5.2E-01 4.3E-02 2.5E-02 ZNF655 LoF 8 1.9E-05 0 4 3.9E-02 3.0E-02 10 5.0E-06 4.6E-04 2.0E-04 RPS16 LoF 3 1.6E-04 2 2 7.4E-01 4.2E-01 3 4.1E-02 7.9E-03 6.4E-03 Table shows genes with P-value < 5.7 × 10−5(High-Mod), 2.2 × 10−4(High), or 2.8 × 10−4(LoF) in the total discovery sample. Results surpassing discovery stage Bonferroni corrected significance thresholds -- P=2.8 × 10−6(High-Mod), 1.1 × 10−5(High), and 1.4 × 10−5(LoF) –are indicated in bold. alocated in APOE region btype of functional variants included in gene-based test Whole exome sequencing study identifies novel rare and common Alzheimer’s-Associated variants. . . ADSP whole-genome sequenced (WGS) family-based dataset. A rare missense variant (rs61756195, MAF =0.001) in STAG3 segregated with disease in three CH families and trended toward association in the case-control study Fig. 2 Manhattan plots showing exome-wide association results for gene-based tests of rare functional variants. The plots show the genebased p-values from the Discovery meta-analysis against their genomic position for association with AD. Each point represents a p-value from SKAT-O test aggregating rare variants (MAF < 5%), by gene, on the basis of predicted functional impact. Only genes with a cumulative minor allele count of ≥10 were included; the minimum p-value from the three adjustment models for either the meta-analysis, European Ancestry (EA), or Caribbean Hispanic (CH) is plotted for each variant. Genes are indicated above points that surpassed our significance threshold for follow-up in tests aggregating only (a) moderate or high impact variants, (b) high impact variants; (c) loss-of-function variants. In each plot, the dotted line indicates the threshold for follow-up: (a) p<5.5 × 10−5,(b) p<6.3 × 10−5,(c) p<2.8 × 10−4, each corresponding to 1 / # genes tested. The dashed line indicates the threshold for exome-wide significance: (a) p<2.7 × 10−6,(b) p<3.1 × 10−6,(c)p< 1.4 × 10−5, each corresponding to 0.05 / # genes tested J. C. Bis et al. (p =0.052) (Supplementary Table S8). Gene-based testing identified a nominal association for GAS2L2 (p =0.049) in EA families (Supplementary Table S9). Rare variants in established genes from GWAS We interrogated our individual variant and gene-based aggregate association tests for 299 previously associated AD genes. Among the SNVs and indels, a total of 1,172 variants with MAF < 0.05 and annotated as either HIGH or MODERATE impact are located within 253 interrogated AD genes (Supplementary Table S10). Five of these variants were at least suggestively significant (p<8.9 × 10−4) in single variant testing. The most significant associations included the TREM2 R47H missense mutation (p =4.8 × 10−12) and ABCA7 frameshift mutation E709fs (p =4.3 × 10−6) which was previously associated with AD in Belgian families [40]. Additional notable signals included variants in SORL1 A528T (p =8.7 × 10−5), which was previously associated with AD in a CH population [15], and ACP2 D353E (p =7.8 × 10−4). Perturbation of murine Acp2 causes lysosomal storage deficits, kyphoscoliosis, cerebellar abnormalities, and ataxia [44,45]. For gene-based tests, we aggregated variants on the basis of annotated function, and examined only genes with more than one contributing variant and a cMAC ≥ 10. Of the 299 AD genes, tests were performed on 281 genes aggregating high or moderate impact variants and86geneslimitedtohighimpactvariants.Among these, 13 unique genes surpassed suggestive significance thresholds for high (p < 1.16 × 10−2)orhigh-moderate (p =3.56 × 10−3) impact variants. The strongest associations were observed for moderate impact variants in TREM2 (p =4.81 × 10−12)andSORL1 (p =8.68 × 10−5), and a high impact variant in ABCA7 (p =4.33 × 10−6). Other noteworthy signals included moderate impact variants in NUP88 (p =4.63 × 10−4)andACP2 (p = 7.80 × 10−4). Discussion Our WES study, the largest for AD conducted to date, identified novel associations with variants in three genes not previously implicated in AD including one common nearly exome-wide significant variant each in IGHG3 (p =9.8 × 10−7) and STAG3 (p =8.8 × 10−7), and one rare exome-wide significant variant in AC099552.4 (p =1.2 × 10−7). We also observed a gene-wide significant association with ZNF655 in a gene-based test including nine highimpact rare variants (p =5.0 × 10−6). These results remained significant after multiple test correction and were confirmed in or strengthened by a replication sample comprised of four independent datasets, with the exception of the variant in AC099552.4 which was invariant in the replication samples. We also confirmed associations with common and rare variants in several previously established AD genes including ABCA7, APOE, HLA-DPA1, MS4A6A, PILRA, SORL1 and TREM2. ZNF655 is expressed in brain and encodes the Vavinteracting Krüppel-like factor 1 [46]. Krüppel-like factors (KLFs) are zinc finger-containing transcription factors that regulate diverse biological processes, including proliferation, differentiation, growth, development, survival, and responses to external stress [47]. Several KLFs have been shown to participate in neuronal morphogenesis and to control the regenerative capacity of neurons in the central nervous system. AC099552.4 is a long non-coding RNA, an abundant class of RNA sequences which regulate gene transcription and expression [48] and impact neuronal development, neuroplasticity, and cognition [49]. Noncoding RNA-dependent regulation affecting AD-related processes has been demonstrated for SORL1 [50] and in a triple transgenic model of AD [51]. IGHG3 encodes immunoglobulin heavy constant gamma 3 and is a member of the IgG family for which antibodies have been shown to cross-react with fibril and oligomer amyloid-βaggregates [52] leading to speculation that Immunoglobulin GM (γmarker) genes contain functional risk and protective factors for AD [53]. The antiamyloidogenic activity of IgG appears to be an inherent property of free human IgG heavy chains [54]. Recent analysis of structural variants in whole genome sequence data for 578 members of 101 families with multiple AD subjects included in the ADSP [26] yielded additional evidence supporting IGHG3 as an AD risk locus. A total of nine distinct deletions in the IGH region were identified as disproportionately represented in AD cases compared to controls. One of these is a 188 bp deletion that was observed in 35 AD cases and 8 controls and is located 592 bp from the AD-associated SNV (rs12890621) in this study. This deletion eliminates a large portion of IGHG3 intron 2 and exon 3 (reference transcript ENST00000390551), and is Fig. 3 Distribution of high impact, LoF and high-confidence LoF variants grouped by predicted consequence Whole exome sequencing study identifies novel rare and common Alzheimer’s-Associated variants. . . and N-ethylmaleimide-sensitive fusion protein in schizophrenic patients. Neurosci Lett. 2001;305:185–8. 61. Reitz C, Cheng R, Rogaeva E, Lee JH, Tokuhiro S, Zou F, et al. Meta-analysis of the association between variants in SORL1 and Alzheimer disease. Arch Neurol. 2011;68:99–106. 62. Rogaeva E, Meng Y, Lee JH, Gu Y, Kawarai T, Zou F, et al. The neuronal sortilin-related receptor SORL1 is genetically associated with Alzheimer disease. Nat Genet. 2007;39:168–77. 63. Willnow TE, Andersen OM. Sorting receptor SORLA-a trafficking path to avoid Alzheimer disease. J Cell Sci. 2013;126: 2751–60. 64. Gavioli EC, de Medeiros IU, Monteiro MC, Calo G, Romao PR. Nociceptin/orphanin FQ-NOP receptor system in inflammatory and immune-mediated diseases. Vitam Horm. 2015;97:241–66. 65. Abdel-Mouttalib O. Nociceptin/orphanin-FQ modulation of learning and memory. Vitam Horm. 2015;97:323–45. 66. Hashimoto R, Noguchi H, Hori H, Nakabayashi T, Suzuki T, Iwata N, et al. A genetic variation in the dysbindin gene (DTNBP1) is associated with memory performance in healthy controls. World J Biol Psychiatry. 2010;11:431–8. 67. Hashimoto R, Noguchi H, Hori H, Ohi K, Yasuda Y, Takeda M, et al. Association between the dysbindin gene (DTNBP1) and cognitive functions in Japanese subjects. Psychiatry Clin Neurosci. 2009;63:550–6. 68. Hemmi H, Idoyaga J, Suda K, Suda N, Kennedy K, Noda M, et al. A new triggering receptor expressed on myeloid cells (Trem) family member, Trem-like 4, binds to dead cells and is a DNAX activation protein 12-linked marker for subsets of mouse macrophages and dendritic cells. J Immunol. 2009;182:1278–86. 69. Mentrup T, Fluhrer R, Schroder B. Latest emerging functions of SPP/SPPL intramembrane proteases. Eur J Cell Biol. 2017;96: 372–82. 70. Verkerk AJ, Schot R, Dumee B, Schellekens K, Swagemakers S, Bertoli-Avella AM, et al. Mutation in the AP4M1 gene provides a model for neuroaxonal injury in cerebral palsy. Am J Hum Genet. 2009;85:40–52. 71. Toh WH, Tan JZ, Zulkefli KL, Houghton FJ, Gleeson PA. Amyloid precursor protein traffics from the Golgi directly to early endosomes in an Arl5band AP4-dependent pathway. Traffic. 2017;18:159–75. Affiliations Joshua C. Bis1●Xueqiu Jian2●Brian W. Kunkle3●Yuning Chen 4●Kara L. Hamilton-Nelson3●William S. Bush5● William J. Salerno6●Daniel Lancour 7●Yiyi Ma7●Alan E. Renton 8●Edoardo Marcora8,9 ●John J. Farrell7● Yi Zhao10 ●Liming Qu10 ●Shahzad Ahmad11 ●Najaf Amin12 ●Philippe Amouyel 12,13,14 ●Gary W. Beecham3● Jennifer E. Below15 ●Dominique Campion16,17 ●Laura Cantwell10 ●Camille Charbonnier16 ●Jaeyoon Chung 7● Paul K. Crane1●Carlos Cruchaga18 ●L. Adrienne Cupples4,19 ●Jean-François Dartigues20 ●Stéphanie Debette20,21 ● Jean-François Deleuze22 ●Lucinda Fulton23 ●Stacey B. Gabriel24 ●Emmanuelle Genin25 ●Richard A. Gibbs6● Alison Goate 8,9 ●Benjamin Grenier-Boley12 ●Namrata Gupta24 ●Jonathan L. Haines5●Aki S. Havulinna26,27 ● Seppo Helisalmi28 ●Mikko Hiltunen29 ●Daniel P. Howrigan30,31 ●M. Arfan Ikram 11 ●Jaakko Kaprio 26 ● Jan Konrad18 ●Amanda Kuzma10 ●Eric S. Lander24 ●Mark Lathrop32 ●Terho Lehtimäki33 ●Honghuang Lin 34 ● Kari Mattila33 ●Richard Mayeux35 ●Donna M. Muzny6●Waleed Nasser6●Benjamin Neale 30,31 ●Kwangsik Nho36 ● Gaël Nicolas 16 ●Devanshi Patel7●Margaret A. Pericak-Vance3●Markus Perola26,27,37 ●Bruce M. Psaty1,38,39,40 ● Olivier Quenez16 ●Farid Rajabli3●Richard Redon41 ●Christiane Reitz35 ●Anne M. Remes28,42 ●Veikko Salomaa27 ● Chloe Sarnowski4●Helena Schmidt43 ●Michael Schmidt3●Reinhold Schmidt43 ●Hilkka Soininen28,44 ● Timothy A. Thornton45 ●Giuseppe Tosto35 ●Christophe Tzourio20 ●Sven J. van der Lee 11 ●Cornelia M. van Duijn11 ● Otto Valladares10 ●Badri Vardarajan35 ●Li-San Wang10 ●Weixin Wang10 ●Ellen Wijsman46,47 ●Richard K. Wilson23 ● Daniela Witten45,47 ●Kim C. Worley6●Xiaoling Zhang4,7 Alzheimer’s Disease Sequencing Project ● Celine Bellenguez 12 ●Jean-Charles Lambert 12 ●Mitja I. Kurki26,30,31 ●Aarno Palotie26,30,31 ●Mark Daly 24,26,31 ● Eric Boerwinkle6,48 ●Kathryn L. Lunetta 4●Anita L. Destefano4,49 ●Josée Dupuis4●Eden R. Martin3● Gerard D. Schellenberg10 ●Sudha Seshadri19,49,50 ●Adam C. Naj 10 ●Myriam Fornage2,48 ●Lindsay A. Farrer4,7,49,51,52 1Department of Medicine (General Internal Medicine), University of Washington, Seattle, WA, USA 2Institute of Molecular Medicine, McGovern Medical School, University of Texas Health Science Center at Houston, Houston, TX, USA 3John P. Hussman Institute for Human Genomics, Miller School of Medicine, University of Miami, Miami, FL, USA 4Departments of Biostatistics, Boston University School of Public Health, Boston, MA, USA 5Case Western Reserve University, Cleveland Heights, OH, USA 6Human Genome Sequencing Center and Department of Molecular and Human Genetics, Baylor College of Medicine, Houston, TX, USA 7Department of Medicine (Biomedical Genetics), Boston University School of Medicine, Boston, MA, USA 8Department of Neuroscience and Ronald M Loeb Center for Alzheimer’s Disease, Icahn School of Medicine at Mount Sinai, New York, NY, USA 9Department of Genetics and Genomics Sciences, Icahn School of Medicine at Mount Sinai, New York, NY, USA 10 University of Pennsylvania Perelman School of Medicine, Philadelphia, PA, USA J. C. Bis et al. 11 Erasmus University Medical Center, Rotterdam, Netherlands 12 Inserm, U1167, RID-AGE-Risk Factors and Molecular Determinants of Aging-Related Diseases, Lille, France 13 Institut Pasteur de Lille, Lille, France 14 University Lille, U1167-Excellence Laboratory LabEx DISTALZ, Lille, France 15 Department of Medical Genetics, Vanderbilt University Medical Center, Nashville, TN, USA 16 Department of Genetics and CNR-MAJ, Normandie Université, UNIROUEN, Inserm U1245 and Rouen University Hospital, F 76000, Normandy Centre for Genomic and Personalized Medicine, Rouen, France 17 Department of Research, Centre Hospitalier du Rouvray, Sotteville-lès-, Rouen, France 18 Department of Psychiatry, Washington University, St. Louis, MO, USA 19 National Heart, Lung, and Blood Institute’s Framingham Heart Study, Framingham, MA, USA 20 University of Bordeaux, Inserm, Bordeaux Population Health Research Center, team VINTAGE, UMR 1219, F-33000 Bordeaux, France 21 Department of Neurology and Institute for Neurodegenerative Diseases, Bordeaux University Hospital, Memory Clinic, F-33000 Bordeaux, France 22 Centre National de Recherche en Génomique Humaine, Institut François Jacob, Direction de le Recherche Fondamentale, CEA, Evry, France 23 McDonnell Genome Institute, Washington University, St. Louis, MO, USA 24 Broad Institute of MIT and Harvard, Cambridge, MA, USA 25 Inserm UMR-1078, CHRU Brest, Université Brest, Brest, France 26 Institute for Molecular Medicine Finland (FIMM), University of Helsinki, Helsinki, Finland 27 National Institute for Health and Welfare, Helsinki, Finland 28 Institute of Clinical Medicine - Neurology and Department of Neurology, University of Eastern Finland, Kuopio, Finland 29 Institute of Biomedicine, University of Eastern Finland, Kuopio, Finland 30 Program in Medical and Population Genetics and Genetic Analysis Platform, Stanley Center for Psychiatric Research, Broad Institute of MIT and Harvard, Cambridge, MA, USA 31 Psychiatric & Neurodevelopmental Genetics Unit, Massachusetts General Hospital, Boston, MA, USA 32 McGill University and Génome Québec Innovation Centre, Montréal, Canada 33 Department of Clinical Chemistry, Fimlab Laboratories and Finnish Cardiovascular Research Center-Tampere, Faculty of Medicine and Life Sciences, University of Tampere, Tampere, Finland 34 Department of Medicine (Computational Biomedicine), Boston University School of Medicine, Boston, MA, USA 35 Columbia University, New York, NY, USA 36 Indiana University School of Medicine, Indianapolis, IN, USA 37 University of Tartu, Estonian Genome Center, Tartu, Estonia 38 Department of Epidemiology, University of Washington, Seattle, WA, USA 39 Department of Health Services, University of Washington, Seattle, WA, USA 40 Kaiser Permanente Washington Health Research Institute, Seattle, WA, USA 41 Inserm, CNRS, Univ. Nantes, CHU Nantes, l’institut du thorax, Nantes, France 42 Unit of Clinical Neuroscience, Neurology, University of Oulu and Medical Research Center, Oulu University Hospital, Oulu, Finland 43 Department of Neurology, Clinical Division of Neurogeriatrics, Medical University of Graz, Graz, Austria 44 Department of Neurology, Kuopio University Hospital, Kuopio, Finland 45 Department of Statistics, University of Washington, Seattle, WA, USA 46 Department of Medicine (Medical Genetics), University of Washington, Seattle, WA, USA 47 Department of Biostatistics, University of Washington, Seattle, WA, USA 48 School of Public Health, University of Texas Health Science Center at Houston, Houston, TX, USA 49 Departments of Neurology, Boston University School of Medicine, Boston, MA, USA 50 Glenn Biggs Institute for Alzheimer’s and Neurodegenerative Diseases, University of Texas Health Sciences Center, San Antonio, TX, USA 51 Department of Epidemiology, Boston University School of Public Health, Boston, MA, USA 52 Department of Ophthalmology, Boston University School of Medicine, Boston, MA, USA Whole exome sequencing study identifies novel rare and common Alzheimer’s-Associated variants. . .