Rare Variants in PLD3 Do Not Affect Risk for Early-Onset Alzheimer Disease in a European Consortium Cohort.

Rare variants in the phospholipase D3 gene (PLD3) were associated with increased risk for late-onset Alzheimer disease (LOAD). We identified a missense mutation in PLD3 in whole-genome sequence data of a patient with autopsy confirmed Alzheimer disease (AD) and onset age of 50 years. Subsequently, we sequenced PLD3 in a Belgian early-onset Alzheimer disease (EOAD) patient (N = 261) and control (N = 319) cohort, as well as in European EOAD patients (N = 946) and control individuals (N = 1,209) ascertained in different European countries. Overall, we identified 22 rare variants with a minor allele frequency <1%, 20 missense and two splicing mutations. Burden analysis did not provide significant evidence for an enrichment of rare PLD3 variants in EOAD patients in any of the patient/control cohorts. Also, meta-analysis of the PLD3 data, including a published dataset of a German EOAD cohort, was not significant (P = 0.43; OR = 1.53, 95% CI 0.60-3.31). Consequently, our data do not support a role for PLD3 rare variants in the genetic etiology of EOAD in European EOAD patients. Our data corroborate the negative replication data obtained in LOAD studies and therefore a genetic role of PLD3 in AD remains to be demonstrated.

Introduction

Whole‐exome sequencing (WES) in familial late‐onset Alzheimer disease (LOAD) patients identified phospholipase D3 (PLD3; MIM #615698) as a new risk gene, with risk variants doubling risk to develop LOAD [Cruchaga et al., 2014]. Typically, phospholipase D (PLD) hydrolyzes the phosphodiester bond of phosphatidylcholine to produce phosphatidic acid (PA) and free choline. PA is an intracellular lipid mediator of many biological functions attributed to PLD, including calcium mobilization, secretion, vesicle trafficking, glucose transport, mitogenesis, and apoptosis [Jang et al., 2008]. A functional link between the PLDs and AD, particularly PLD1, has been previously established. PLD1 was shown to interact with presenilin 1 (PSEN1), to antagonize the PSEN1‐mediated cleavage of βAPP [Cai et al., 2006a] and to affect its intracellular trafficking [Cai et al., 2006b]. For PLD3, the canonical activity had not yet been identified [Munck et al., 2005], but Cruchaga et al. (2014) demonstrated an inverse correlation between PLD3 and amyloid‐β in vitro, proposing a role for PLD3 in APP processing. The PLD3 variant, p.A442, was associated with reduced PLD3 mRNA expression, corroborating the proposed role of PLD3 in LOAD [Cruchaga et al., 2014]. Another variant, p.V232M, was associated with both disease risk and earlier onset age [Cruchaga et al., 2014], suggesting a possible role for PLD3 variants in risk for early‐onset AD (EOAD).

We analyzed available whole‐genome‐sequencing (WGS) data of 20 AD patients with early onset of disease and WES data of five probands diagnosed with autosomal dominant inherited AD, for genetic variation in known causal and risk genes reported for AD. We identified one missense mutation in PLD3 in a familial EOAD patient with autopsy confirmed AD and onset of disease at the age of 50 years. This finding prompted us to screen PLD3 in cohorts of Belgian EOAD patients and matched control individuals. No evidence for a genetic risk association was obtained. Next, we replicated our PLD3 findings in four European EOAD patient/control cohorts ascertained within the European Early‐Onset Dementia consortium [van der Zee et al., 2013; van der Zee et al., 2014; Cuyvers et al., 2015]. While our studies were in progress, four publications provided evidence that the initial PLD3 risk association, obtained in a family‐based study, could not be replicated in independently ascertained LOAD patient/control cohorts [Lambert et al., 2015; van der Lee et al., 2015; Heilmann et al., 2015; Hooli et al., 2015], questioning a role for rare PLD3 variants in AD risk.

Methods

Belgian Patient/Control Cohorts

The EOAD patient group consisted of unrelated AD patients (n = 286), aged ≤70 years at inclusion, of which n = 123 (43%) were diagnosed at age ≤65 years. The EOAD patients were selected from a larger AD patient group sampled within a prospective study of Belgian dementia patients at the Neurology Department and Memory Clinic of the Hospital Network Antwerp (ZNA) Middelheim and Hoge Beuken, (P.P.D.D. and S.E.) [Engelborghs et al., 2003; Engelborghs et al., 2006] and the Department of Neurology and the Memory Clinic of the University Hospitals of Leuven, Leuven, (R.V. and M.V.) [Bettens et al., 2010], as well as through the neurology centers of the clinical partners within the Belgian Neurology (BELNEU) consortium. In the EOAD patient cohort, average onset age was 63.6 ± 6.3 years (range 37–70; in addition, two autosomal dominant AD probands had a later onset age of 76 and 77 years, respectively), 58% were women, and 37 of the patients received a definite diagnosis of AD at autopsy. Consensus diagnosis of possible, probable or definite AD was given by at least two neurologists based on the National Institute of Neurological and Communicative Disorders and Stroke (NINCDS) and the Alzheimer's Disease and Related Disorders Association (ADRDA) or the National Institute on Aging‐Alzheimer's Association (NIA‐AA) diagnostic criteria [McKhann et al., 1984; McKhann et al., 2011; Hyman et al., 2012]. Each AD patient underwent a neuropsychological examination, including Mini‐Mental State Examination (MMSE) [Folstein et al., 1975] and structural imaging, and a subset of patients was also examined by functional neuroimaging [Bettens et al., 2010].

The Flanders‐Belgian control cohort consisted of 319 unrelated older individuals with a mean age at inclusion of 68.7 ± 12.2 years and 64% women. Subjective memory complaints and neurological or psychiatric antecedents as well as a familial history of neurodegeneration were ruled out by means of an interview. Cognitive screening was performed using the MMSE (inclusion cutoff score ≥26) [Folstein et al., 1975] or the Montreal Cognitive Assessment (inclusion cutoff score ≥25), a 30‐point test known to be sensitive for Mild Cognitive Impairment (MCI) and early AD [Nasreddine et al., 2005]. The majority (n = 233) of the control individuals consisted of community‐dwelling volunteers and spouses of patients, who were recruited in the hospital upon neurological examination of their affected spouse.

All participants and/or their legal guardian in case of dementia patients, signed a written informed consent form for participation in clinical and genetic studies. The clinical study protocols and the informed consent forms for patients and control individuals were approved by the Ethics Committee of the respective hospitals and/or universities involved in the sampling in Belgium through partners of the BELNEU consortium. The genetic study protocol was approved by the Ethics Committee of the University of Antwerp and the University Hospital of Antwerp, Belgium.

EU Patient/Control Cohorts

EOAD patients and control individuals were ascertained by EU partners in the EU EOD consortium [van der Zee et al., 2013; van der Zee et al., 2014; Cuyvers et al., 2015]. DNA, medical, and demographic information were available of 946 EOAD patients with onset age of ≤65 years. For this study, the individual countries contributed 313 (Spain), 108 (Portugal), 194 (Italy), 172 (Sweden), 98 (Germany), 54 (Greece), and seven (Czech Republic) EOAD patients. The German patient cohort overlaps in part with the EOAD cohort screened for PLD3 variations by Schulte and colleagues [Schulte et al., 2015]. The EOAD patients were diagnosed following the NINCDS‐ADRDA International Working Group criteria [McKhann et al., 1984; McKhann et al., 2011]. Pathological diagnosis of definite AD patients was based on current diagnostic criteria [Montine et al., 2012]. A definite AD diagnosis was available for 40 patients from Spain, two patients from Portugal, five patients from Italy, two patients from Sweden, and seven patients from Czech Republic.

The total EU control group consisted of 1,209 persons contributed by partnering countries Spain (303), Portugal (126), Italy (428), Sweden (345), and Czech Republic (7). Control persons had normal cognition based on MMSE (inclusion cutoff score ≥26). All EU participants signed informed consent forms that were approved by the ethics committee of the local hospitals or sampling sites.

PLD3 Variation in WGS and WES data of Belgian EOAD Patients

For WGS studies, we had selected 20 patients from the EOAD cohort with an early‐onset age (mean 56.7 years ± 6.6 years; range 37–65 years), of whom 15 were negative for the ε4 allele of the apolipoprotein E gene (APOE). An additional five probands of autosomal dominant families (mean 61.4 years; range 50–77 years) were selected for WES. For all 25 patients, pathological mutations in APP, PSEN1, and PSEN2 were excluded.

WGS was outsourced to Complete Genomics (Mountain View, CA) making use of unchained combinatorial probe anchor ligation chemistry on self‐assembling DNA nanoballs [Drmanac et al., 2010]. Raw data reads were assembled to the reference genome (National Center for Biotechnology Information build 36 or Genome Reference Consortium Human Build 37/hg19 (GRCh37/hg19)). Sequencing alignment and variant calling were performed at Complete Genomics (Complete Genomics, whereas data annotation and analysis were performed by us with the GenomeComb package [Reumers et al., 2012]. WES data were generated with various enrichment kits and sequencing was performed on the Illumina® NextSeq (http://www.illumina.com) or on the Applied Biosystems SOLiD 5500 (http://www.solid.appliedbiosystems.com) systems. Sequence manipulation and variants calling from the Illumina® NextSeq data were performed with GenomeComb [Reumers et al., 2012], which integrates the pipeline for NGS data analysis. Sequencing data obtained from the SOLiD 5500 system where analyzed using LifeScope™ Genomic Analysis Software (http://www.lifetechnologies.com/lifescope), annotated and analyzed using the GenomeComb package [Reumers et al., 2012].

Sequence data of PLD3 (RefSeq NM_001031696.3) coding exons (including intron‐exon boundaries) were extracted from the annotated WGS and WES data. For all sequenced samples, at least 95% of the coding DNA sequence (CDS) of PLD3 was covered ≥20×. Coding and splice variants with a coverage ≥20× and a minor allele frequency (MAF) lower than 1% in the Exome Variant Server (EVS; http://evs.gs.washington.edu/EVS/) were selected.

PLD3 Sequencing in EOAD Patient/Control Cohorts

In the Belgian cohort, the CDS of PLD3 (RefSeq NM_001031696.3) was analyzed by Sanger sequencing (BigDye Terminator Cycle Sequencing kit v3.1 on the ABI 3730 DNA Analyzer (both Applied Biosystems®, Waltham, Massachusetts, USA). Sequences were analyzed for genetic variations using SeqManII (DNASTAR, Madison, WI) and/or novoSNP [Weckx et al., 2005] software packages. All DNAs included in the study were successfully screened for at least 85% of the CDS of PLD3.

To screen the EU EOAD cohorts, we designed an amplicon‐based gene assay for massive parallel sequencing of all coding exons of PLD3. PLD3 exonic regions were taken from Gencode human gene annotation set release v19 [Harrow et al., 2012]. A minimum of 15 nt intronic sequence flanking each exon was included to cover splice sites. Target regions larger than 500 nt in size were divided into contiguous segments of less than 500 nt to anticipate sequence read coverage in subsequent massive parallel sequencing experiments. Adjacent target regions located within a 500 bp segment were combined to form a single target region. PCR Primers flanking each target were designed using mPCR (Multiplicom N.V., Niel, Belgium), a software tool to design multiplex PCR panels for specific and uniform target amplification, built around the Primer3 program [Koressaar and Remm, 2007]. Together, all targets were covered by 19 amplicons in three multiplex PCR reactions. Target‐specific primer sequences were extended with tag sequences (5′‐TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG‐Fwd and 5′‐GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAG‐Rev) for subsequent indexing and sequencing. First, optimal annealing temperature and relative amounts of individual PCR primers were optimized for uniform amplification of each target in the multiplex reaction. Then, multiplex PCR reactions were carried out on 20 ng genomiphied DNA (Illustra GenomiPhi V2; Thermo Fisher Scientific Inc., Waltham, Massachusetts, USA). The amplification quality and efficiency of each reaction was verified by fragment analysis on an ABI 3730 automated sequencer, and the amplification products of all multiplex PCR reactions of each individual were pooled aiming at equimolar representation of the amplicons in the mixture. Per‐individual indexing was performed using the Nextera indexing system (Illumina®) allowing the pooling of up to 1,536 distinct index combinations. After indexing, amplicons were pooled in three separate sequencing libraries. Massive parallel sequencing was performed on a MiSeq sequencer (Illumina®) using V2 chemistry to generate 2 × 250 paired‐end sequencing reads. After sample demultiplexing, sequence reads were mapped using the Burrows‐Wheeler Aligner [Li and Durbin, 2009; Li and Durbin, 2010] to a minigenome consisting of the combined amplicon sequences extracted from the human genome reference sequence hg19. Sequence variants were called using the Genome Analysis Toolkit [McKenna et al., 2010; DePristo et al., 2011]and annotated using the GenomeComb variant annotation pipeline [Reumers et al., 2012].

Variants with a read depth below 15×, with Genoqual value below 99 or with an imbalanced wild‐type/variant read depth (cutoff > 4) were considered false calls. All remaining rare (MAF < 1%) variants, predicted to affect protein sequence and splicing variants, were validated on genomic DNA by Sanger sequencing. We also selected the silent p.A442 variant because of the previous reported AD risk association [Cruchaga et al., 2014], and included it in downstream statistical analysis.

The variants reported are based on the RefSeq NM_001031696.3. The nucleotide positions refer to the cDNA sequence and the nucleotide numbering uses +1 as the ATG translation initiation codon in the reference sequence, with the initiation codon as codon 1. The variants identified in this study have been submitted to the Leiden Open Variation Database (LOVD) Locus Specific Database (LSDB) (http://databases.lovd.nl/shared/genes/PLD3).

Statistical Analyses

A rare variant (MAF < 1%) burden analysis was performed by collapsing of rare variant alleles and comparing the overall frequency of rare variant alleles between patients and control individuals using χ 2 statistics. Odds ratios (OR) and 95% confidence intervals were computed. A two‐sided P value <0.05 was considered significant. Fixed effect (Mantel–Haenszel) meta‐analysis was performed based on the allelic data of the Belgian and EU EOAD cohorts. Because of too small sample size or lack of a matched control group, we did not include the Czech, German, and Greek cohorts in the association analysis. Meta‐analysis was performed on the PLD3 data obtained in the Belgian, Spanish, Portuguese, Italian, and Swedish patient/control cohorts plus the one published German study of PLD3 in EOAD [Schulte et al., 2015]. Mantel–Haenszel summary P value, OR, Higgins I 2, and Cochran's Q test for heterogeneity were computed in R using the library epiR, version 0.9–62.

Characterization of PLD3 Splicing Variants

Bioinformatic prediction analysis of splicing variants was performed using both FSPLICE 1.0 (http://www.softberry.com/; Softberry, Inc., Mount Kisco, NY, USA) and NNSPLICE 0.9 (http://www.fruitfly.org/seq_tools/splice.html) [Reese et al., 1997]. Both programs search for splicing sites in the input sequence. The analysis was performed on both wild‐type and mutated sequences using both tools and the comparison between the outputs was performed to identify differences in the splicing sites identification (Supporting Information).

Total RNA was extracted from lymphoblast cell lines using the RiboPureTM kit, Ambion® (Life TechnologiesTM, Waltham, Massachusetts, USA) followed by DNase treatment, performed using Ambion® TURBO DNA‐freeTM kit (Life TechnologiesTM). Complementary DNA was prepared from the total RNA, using SuperScript® III First‐Strand Synthesis System for RT‐PCR, InvitrogenTM (Life TechnologiesTM) with random hexamer primers including a control reaction without retro‐transcriptase enzyme. cDNA amplification was carried out using Titanium® Taq DNA Polymerase (Clontech Laboratories, Inc., Mountain View, CA, USA). PCR product agarose gel band extraction was performed using QIAquick Gel Extraction Kit (http://www.qiagen.com) following the manufacturer guidelines. The purified amplification products were sequenced (BigDye Terminator Cycle Sequencing kit v3.1; Applied Biosystems®) or cloned into pCRTM2.1‐TOPO® vector following manufacturer instructions and the purified vector was sequenced and analyzed on an ABI 3730 DNA Analyzer (Applied Biosystems®). Sequence analysis was performed using SeqManII (DNASTAR).

Results

PLD3 Variants in the Belgian EOAD Patient/Control Cohort

Analysis of the WGS data of 20 unrelated EOAD patients and WES data of five probands of autosomal dominant AD families revealed one rare PLD3 variant, p.N284S, in a familial EOAD patient with an onset age of 50 years and an APOE ε3ε3 genotype. Review of the clinical course of the disease in the carrier showed no unusual characteristics. The amyloid β T‐tau and P‐tau biomarker profile in cerebrospinal fluid was compatible with AD, and autopsy after a disease duration of 7 years confirmed the diagnosis of AD (Braak and Braak stage V‐VI, Montine classification A3B3C3). Based on heteroanamnesis, both parents had late‐onset dementia >90 years and a paternal uncle had cognitive impairment at younger age. DNA of family members was not available for genetic testing of PLD3.

The sequencing of PLD3 in the 261 remaining Belgian patients, identified four missense mutations and one splicing variation in eight patients (Fig. 1), resulting in a cumulative rare variant frequency in Belgian EOAD patients of 3.15% (nine out of 286) (Table 1). This included one missense, p.P173S, and one splicing variant, c.1285+1G>A (p.K396Efs*6), that had not yet been reported. In 10 control individuals, we observed seven different non‐synonymous mutations (10 out of 319, 3.13%) of which three were new (Table 1). Of the reported LOAD‐associated variants [Cruchaga et al., 2014], we identified p.V232M in one patient and four control persons (P = 0.2). The p.A442 variant was observed in seven patients (seven out of 286, 2.4%) and 10 control individuals (10 out of 319, 3.13%) (P = 0.6). Rare variant gene‐based burden analysis did not show a significant enrichment of PLD3 rare variant alleles in AD patients compared with control individuals (P = 1.0, relative risk (OR) = 1.00, 95% CI = 0.40–2.49).

Figure 1

Location of PLD3 variants. Schematic representation of rare missense and splicing variants identified in the PLD3 (RefSeq NM_001031696.3), in the Belgian and EU EOAD cohorts. For completeness, the p.A442 variant is also reported. Variations in red are observed in patients‐only, in blue in control persons‐only and in black those identified in both patients and control individuals (Table 1). Variants marked with an asterisk (*) were previously reported by Cruchaga et al. (2014). The exon numbering follows the RefSeq NM_001031696.3.

Table 1

Overall Rare Variants Identified in PLD3

				Belgium		Spain		Portugal		Sweden		Italy		Germany		Czech republic		Greece
Total alleles (AD and C)				AD n = 572, C n = 638		AD n = 626, C n = 606		AD n = 216, C n = 252		AD n = 344, C n = 690		AD = 388, C n = 856		AD n = 196, C n.a.		AD n = 14, C n = 14		AD = 108, C n.a.
dbSNP	cDNA	Protein		Alleles	MAF	Alleles	MAF	Alleles	MAF	Alleles	MAF	Alleles	MAF	Alleles	MAF	Alleles	MAF	Alleles	MAF
rs142070038	c.187G>A	p.G63S	AD	0		2	0.32	2	0.93	0		2	0.52	0		0		2	1.85
			C	1	0.16	0		0		0		1	0.12	n.a.		0		n.a.
rs138674695	c.226C>G	p.P76A	AD	0		0		0		0		1	0.26	0		0		0
			C	1	0.16	0		0		2	0.29	0		n.a.		0		n.a
rs374184677	c.475G>A	p.V159M	AD	0		1	0.16	0		0		1	0.26	0		0		0
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.484C>T	p.R162C	AD	0		0		0		2	0.58	0		0		0		0
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.517C>T	p.P173S	AD	1	0.17	0		0		0		0		0		0		0
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.524C>G	p.A175G	AD	0		0		1	0.46	0		0		0		0		0
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.562C>T	p.R188C	AD	0		0		1	0.46	0		0		0		0		0
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.665G>A	p.R222H	AD	0		0		0		0		0		0		0		0
			C	0		0		1	0.40	0		0		n.a.		0		n.a.
rs145999145	c.694G>A	p.V232M	AD	1	0.17	1	0.16	0		2	0.58	0		2	1.02	0		0
			C	4	0.63	1	0.17	0		1	0.14	4	0.47	n.a.		0		n.a.
n.a.	c.725G>A	p.R242Q	AD	0		0		0		0		0		0		0		0
			C	1	0.16	0		0		0		0		n.a.		0		n.a.
n.a.	c.746A>G	p.E249G	AD	0		0		0		0		0		0		0		0
			C	0		0		0		0		1	0.12	n.a.		0		n.a.
rs144312764	c.814C>T	p.R272C	AD	1	0.17	0		0		0		0		0		0		0
			C	1	0.16	0		0		0		0		n.a.		0		n.a.
rs200274020	c.851A>G	p.N284S	AD	1	0.17	2	0.32	0		0		0		1	0.51	0		0
			C	0		3	0.50	0		0		0		n.a.		0		n.a.
rs368737000	c.878C>T	p.A293V	AD	0		0		0		0		0		1	0.51	0		0
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.890C>T	p.P297L	AD	0		0		0		0		0		0		0		0
			C	1	0.16	0		0		0		0		n.a.		0		n.a
rs146083475	c.899G>A	p.C300Y	AD	3	0.52	3	0.48	0		0		1	0.26	0		0		1	0.93
			C	0		0		0		0		3	0.35	n.a.		0		n.a.
n.a.	c.923T>C	p.L308P	AD	0		0		0		0		0		0		0		0
			C	1	0.16	0		0		0		0		n.a.		0		n.a.
n.a.	c.1019+2T>A (out of frame ex10 del)	p.S294Vfs*6	AD	0		0		0		0		0		0		0		0
			C	0		0		0		1	0.14	0		n.a.		0		n.a.
rs370488565	c.1072G>A	p.V358I	AD	1	0.17	0		0		0		0		0		0		0
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.1276A>G	p.T426A	AD	0		0		0		0		0		0		0		2	1.85
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.1285+1G>A (out of frame ex12del)	p.K396Efs*20	AD	1	0.17	0		0		0		0		0		0		0
			C	0		0		0		0		0		n.a.		0		n.a.
n.a.	c.1285G>A	p.G429R	AD	0		0		0		0		0		0		0		0
			C	0		0		0		1	0.14	0		n.a.		0		n.a
rs4819	c.1326G>A	p.A442	AD	7	1.22	11	1.76	2	0.93	3	0.87	1	0.26	5	2.55	1		1
			C	10	1.57	4	0.66	1	0.40	9	1.30	6	0.70	n.a.		1		n.a.

Missense, splicing and rs4819 p.A442 variant (bottom of table) identified in PLD3 in AD patients and control (C) individuals from Belgium, Spain, Portugal, Sweden, Italy, and Czech Republic and in EOAD patients from Germany and Greece. The variants reported are based on the RefSeq NM_001031696.3. The nucleotide positions refer to the cDNA sequence and the nucleotide numbering uses +1 as the ATG translation initiation codon in the reference sequence, with the initiation codon as codon 1. Minor allele frequency (MAF%) is reported for the single variant identified per cohort. n.a., not available. For the splicing variants, the predicted protein change is shown.

PLD3 Replication Analysis in the EU Cohorts

To replicate our finding and increase statistical power, we investigated PLD3 in 946 EOAD patients and 1,209 age‐ and origin‐matched control individuals ascertained by partners of the EU EOD consortium in Spain, Portugal, Italy, Sweden, Germany, Greece, and Czech Republic. In the overall patient/control cohort, we identified 15 genetic variants with a MAF <0.01, 14 missense and one splicing variant (Fig. 1), in a total of 50 participants of which 31 (62%) patients. Six variants—p.V159M, p.R162C, p.A175G, p.R188C, p.A293V, and p.T426A—were present in patients only, and four variants—p.R222H, p.E249G, c.1019+2T>A (p.S294Vfs*6), and p.G429R, in control persons only (Table 1). The LOAD‐associated p.A442 variant was found in 45 persons of which 24 (53.3%) patients. One control individual from Spain carried both the p.N284S and the p.A442 variant.

Rare variant gene burden analysis did not show an enrichment for rare variants in patients when comparing missense and splicing variants in patients against control individuals in separate cohorts nor when the different patient/control cohorts of Belgium, Spain, Portugal, Italy, and Sweden were pooled in a meta‐analysis including one published German study on EOAD (Schulte et al., 2015) (P = 0.43; summary OR = 1.53, 95% CI 0.60–3.31) (Table 2). The LOAD‐associated p.V232M variant was found in the Spanish patient/control cohort in one patient (one out of 313, 0.32%) and one control (one out of 303, 0.33%) (P = 0.98), in the Swedish patient/control cohort in two patients (two out of 172, 1.16%) and one control (one out of 345, 0.29%) (P = 0.22), in the Italian patient/control cohort in four Italian control individuals (four out of 428, 0.93%) and in two of the German patients (Table 1). Meta‐analysis did not show association of p.V232M with disease status (P = 0.44; OR = 0.79, 95% CI 0.07–3.11). Similarly, the p.A442 variant did not show association across cohorts nor when meta‐analysis was performed (P = 0.89; OR = 1.13, 95% CI 0.40–2.22) (Table 2). Because of the functional effect of the p.A442 variant on splicing [Cruchaga et al., 2014], we included this variant together with the missense and the splicing variants in the gene burden and meta‐analysis. Also in this case, we did not identify association with disease status in the meta‐analysis (P = 0.57; OR = 1.32, 95% CI 0.65–2.20) (Fig. 2). In the Spanish patient/control cohort, however, a gene burden positive association was observed when the p.A442 variant was included (P = 0.03; OR 2.47, 95% CI 1.07–5.64) (Fig. 2). The investigated cohorts did not show evidence of heterogeneity (I 2 range 0%–21.08%; Q P value > 0.27). Random‐effects meta‐analysis, according to DerSimonian and Laird, confirmed absence of association (Supporting Information).

Table 2

Meta‐Analysis of PLD3 Rare Variants

	Total alleles (N)		Missense and splicing variants			p.A442 variant
Country	AD	C	AD (MAF%)	C (MAF%)	OR [95% CI]	AD (MAF%)	C (MAF%)	OR [95% CI]
Belgium	572	638	9 (1.57)	10 (1.57)	1.00 [0.40–2.49]	7 (1.22)	10 (1.57)	0.78 [0.29–2.06]
Spain	626	606	9 (1.44)	4 (0.66)	2.20 [0.67–7.17]	11 (1.76)	4 (0.66)	2.69 [0.85–8.50]
Portugal	216	252	4 (1.85)	1 (0.40)	4.74 [0.53–42.69]	2 (0.93)	1 (0.40)	2.35 [0.21–26.05]
Sweden	344	690	4 (1.16)	5 (0.72)	1.61 [0.43–6.04]	3 (0.87)	9 (1.30)	0.67 [0.18–2.47]
Italy	388	856	5 (1.29)	9 (1.05)	1.23 [0.41–3.70]	1 (0.26)	6 (0.70)	0.37 [0.04–3.05]
Schulte et al. (2015)	278	2,912	3 (1.08)	18 (0.62)	1.75 [0.51–5.99]	5 (1.80)	39 (1.34)	1.35 [0.53–3.45]
Summary	2,424	5,954	34 (1.40)	47 (0.79)	1.53 [0.60–3.31]	29 (1.20)	69 (1.16)	1.13 [0.40–2.22]
P MH value					0.43			0.89
Q P value					0.78			0.36
I 2 (%)					0			8.89

Meta‐analysis was performed on the five patient/control cohorts screened in this study together with the published data on EOAD [Schulte et al., 2015]. Rare (MAF <1%) non‐synonymous and splicing alleles were included in the rare variant gene burden analysis. Total allele counts and frequencies of rare non‐synonymous and splicing variants as well as the p.A442 variant (c.1326G>A, RefSeq NM_001031696.3) identified in AD patients and control (C) individuals are shown. The P value is given for fixed‐effect Mantel–Haenszel test (pMH) and Higgins I 2 (%) and Cochran's Q test (P value) for heterogeneity are shown for the meta‐analysis.

Figure 2

Forest plot for meta‐analysis of rare variant association. Forest plot for meta‐analysis of the rare missense and splicing variants, and the p.A442 variant. Fixed‐effects meta‐analysis is based on raw allelic data. Five different countries (Belgium, Spain, Portugal, Sweden, and Italy) in addition to the German cohort published by Schulte et al. (2015), are included in the statistical analysis. Summary statistics, fixed‐effects P value (P MH‐value) and odds ratio (OR) including 95% confidence interval (CI) is given in the figure. Evidence of heterogeneity is tested using Cochran's Q test (P value) and the Higgins I 2 statistics (%).

PLD3 Splicing Variants Lead to Exon Skipping

Two carriers of a PLD3 splicing variant were identified: c.1019+2T>A (hg19, chromosome 19, g.40880529) in a Swedish control individual aged 61 years at inclusion, and c.1285+1G>A (hg19, chromosome 19, g.40883793) in a Belgian AD patient with onset age 69 years. Two independent software tools predicted that the c.1019+2T>A variant is likely to abolish the canonical splice donor site of exon 10 (alias exon 8 in [Cruchaga et al., 2014]). This will cause skipping of exon 10 and the predicted introduction of a premature stop codon after 298 amino acids (p.S294Vfs*6). Cells were not available to validate whether the mutant transcript is degraded by the nonsense‐medicated mRNA decay (NMD) control system. The c.1285+1G>A variant is located at the first base of the canonical splice donor site of exon 12 (alias exon 10) and is predicted to abolish the donor splice site of exon 12, causing out‐frame skipping of exon 12 and the premature introduction of a stop codon in exon 13, the last exon of PLD3, leading to the predicted protein change p.K396Efs*20. We produced cDNA from lymphoblast cells of the patient carrier and confirmed exon 12 skipping by Sanger sequencing and agarose gel electrophoresis (Fig. 3). In addition to the two expected alleles, we observed a third band that was about 100 bp longer than the normal allele. This longer band was best visible after amplification of exon 11 to 3′UTR PLD3 cDNA fragment even though also detectable when amplifying the full‐length PLD3 cDNA. In the same experiment, it was also faintly visible in control individuals (Fig. 3). Sequencing of the bands extracted from the 4% agarose gel of the c.1285+1G>A carrier confirmed the absence of exon 12 in the shorter fragment (lower band) and the wild‐type sequence in the normal fragment (middle band). The sequence of the longer fragment, which was cloned in a sequencing vector after extraction form the agarose gel, indicated the presence of the intron 12 sequence and the c.1285+1G>A variation as the result of intron retention in the absence of an active exon 12 splice donor site (Fig. 3).

Figure 3

PLD3 exon 12 skipping and NMD escape. Representation of PLD3 exon 12 skipping mechanism. Exon counting is based on RefSeq NM_001031696.3. In Cruchaga et al. (2014), exon 12 corresponds to exon 10. Agarose gel picture (top) shows the results of the exon 11 ‐ 3′UTR fragments amplification, which was carried out using a primer in exon 11 (5′‐ACGAGCGTGGCGTCAAG‐3′) [Cruchaga et al., 2014] in combination with a reverse primer located in the 3′‐UTR (5′‐GGGGTGGAGGTAGAGGTGGG‐3′). The wild‐type amplicon is 563 bp in size and in case of exon 12 skipping (100 bp) the expected amplicon length is 463 bp, in case of intron retention the product is ∼660 bp. The full‐length PLD3 cDNA amplification product is reported (bottom gel). Amplification carried out using forward primer 5′‐ATCCCCCAGCCTTGAGGGAAG‐3′ and reverse primer 5′‐GGGGTGGAGGTAGAGGTGGG ‐3′). The wild‐type PCR fragment is 1,644 bp and the two additional bands are visible. Electropherograms, obtained after sequencing of the agarose gel extraction product, show the normal transcript and exon 12 skipping. The longer product shows retention of intron 12. In the right bottom panel, the splicing mechanism is showed with the production of the three alternative transcripts (a) wild‐type, (b) exon 12 skipping, and (c) longer transcript with intron 12 retention). A blue triangle in the exon 13 of the aberrant transcript indicates the introduction of a premature termination codon (PTC).

Discussion

Rare variants in PLD3 have been associated with an increased risk of developing LOAD in one study [Cruchaga et al., 2014]. Moreover, in that study, the missense variant p.V232M was also associated with onset age and the silent variant, p.A442, was associated with LOAD risk and lower levels of both PLD3 total mRNA and of PLD3 exon 13 (alias exon 11)‐containing transcripts. The p.A442 variant was predicted to affect alternative splicing by changing the consensus sequence of an exonic splicing enhancer binding site [Cruchaga et al., 2014]. PLD3 is highly expressed in hippocampus and cerebral cortex, critical regions for AD pathology, and it was shown that senile plaques in the frontal cortex containing dystrophic neurites have intense PLD3 immunoreactivity [Satoh et al., 2014]. In case of decreased PLD3 protein expression, APP processing is altered and more extensive amyloid‐β pathology is detected [Cruchaga et al., 2014]. We investigated, in a first stage, the WGS and WES data of a selection of high genetic load Belgian EOAD patients, collected and prioritized for their potential of harboring high penetrant, disruptive variants responsible for the onset of AD at younger age. The identification of the p.N284S variant in PLD3 in one familial EOAD patient with an onset age of 50 years and an autopsy diagnosis of AD, prompted us to analyze a Belgian cohort of EOAD patients (N = 261) without a mutation in the three AD genes APP, PSEN1, and PSEN2, and matched control individuals. Although we did identify additional genetic variations in PLD3, both new and reported ones, the rare variant burden analysis was not significant. To replicate this negative finding, we analyzed patient/control cohorts of different European countries ascertained by partners from the European EOD consortium, in total, 946 EOAD patients and 1,209 control individuals. In total, we identified 22 rare variants with a MAF <1%, 20 missense and two splicing variants. Nonetheless, we did not identify an enrichment of rare variants in patients compared with control individuals, neither when we investigated gene burden or single variants that were previously associated with LOAD per country, nor in the meta‐analysis. Only the Spanish cohort showed an association signal but, solely when the rare non‐synonymous variants were pooled with the p.A442 variant, which was not significant in the single variant analysis, suggesting a false positive association in the Spanish cohort. This is supported by the results of the meta‐analysis which did not show any association with rare variants in general in PLD3 and EOAD risk. Amongst the rare variants identified in the current study, two were novel splice site variants. One was identified in an AD patient and one in a control individual. We showed that the splice variant identified in the AD patient does not lead to NMD because the premature stop codon is introduced in the last exon of the transcript. Moreover, Schulte et al. (2015) also reported an exon deletion variant in one control individual, suggesting that this mechanism may not be related to AD pathology.

While the current study was in progress, different attempts to replicate the genetic association of rare variants in PLD3 with LOAD risk in family‐based studies, case–control cohorts, and population cohorts were unsuccessful [Lambert et al., 2015; van der Lee et al., 2015; Heilmann et al., 2015; Hooli et al., 2015]. Our data in EOAD cohorts corroborate these negative findings. Of note, both the p.V232M and the p.A442 variants from the initial study [Cruchaga et al., 2014] were not associated with AD risk in any of the EOAD cohorts we investigated as well as in the other LOAD studies [Lambert et al., 2015; van der Lee et al., 2015; Heilmann et al., 2015; Hooli et al., 2015]. So far, only one study reported a nominal excess of PLD3 rare variants when EOAD patients were pooled with LOAD (P = 0.03) but the gene‐burden association was not significant when EOAD patients were investigated separately (P = 0.54, n = 139) [Schulte et al., 2015]. Notably, also in this study, neither p.V232M nor p.A442 were associated with EOAD or LOAD risk [Schulte et al., 2015]. On the other hand, they identified a silent variant, p.I364, present exclusively in LOAD patients with a MAF of 1.05% in 15 heterozygous and four homozygous carriers that seems to contribute significantly to the observed association [Schulte et al., 2015]. We included the EOAD cohort screened by Schulte et al. (2015) in our meta‐analysis with the European cohorts screened in this study and overall we did not identify a significant association of rare PLD3 variants when patients were compared with control individuals.

The authors of the original PLD3 association finding in LOAD [Cruchaga and Goate, 2015a; Cruchaga and Goate, 2015b], attributed the failure of the replication studies [Lambert et al., 2015; van der Lee et al., 2015; Heilmann et al., 2015] to population stratification and/or regional differences in variant MAFs. The latter factors might have a major influence in replication studies of low‐frequency variants [Cruchaga and Goate, 2015a]. We indeed noticed that the appearance and MAFs of rare variants in the different European cohorts was dissimilar. For example, we did not observe the previously, nominally associated p.M6R variant [Cruchaga et al., 2014] in any of the cohorts we screened. Nevertheless, each patient/control cohort in our study was ethnically matched, arguing against a false‐negative finding due to cohort stratification.

In conclusion, our study investigated the role of rare variants in PLD3 in an extended cohort of European EOAD patient and matched control individuals that were ascertained in the BELNEU and EU EOD consortia. Our study was not limited to the rare variants published [Cruchaga et al., 2014], but included a complete analysis of the coding region to capture all rare variants in patients and matched control persons. The negative replication data we obtained were in line with those reported in four studies of LOAD [Lambert et al., 2015; van der Lee et al., 2015; Heilmann et al., 2015; Hooli et al., 2015]. Our findings in EOAD, reinforce the notion that PLD3 unlikely contributes to genetic risk for AD, indifferent of age of onset of the disease. Moreover, if PLD3 variants were involved in earlier onset age of AD symptoms, as was suggested in the initial study, one could expect that a cohort of EOAD patients would be enriched for these variants compared with LOAD patients. This was clearly not the case in our study making it unlikely that there is a correlation between PLD3 variations and earlier‐onset age. Together, the contradictory data indicate that taken the negative replication data, we should be careful to attribute a biological role to PLD3 in AD pathogenesis.

Supporting Information

Click here for additional data file.