Introduction

Autism spectrum disorders (ASD) have a strong genetic etiology. The last decade of research in psychiatric genetics has led to groundbreaking discoveries, including the identification of ASD-associated sequence and structural variants in a variety of genes implicated in brain development—particularly in synapse formation, maintenance and function. However, the mutations identified so far explain only a fraction of ASD cases, underscoring the need to continue the search for missing genetic links in ASD.

Wnt intercellular signaling refers to a set of key developmental pathways that help to organize the nervous system. Among other roles, Wnt signaling influences subdivision of the anterior neural tube into major brain regions (such as forebrain, midbrain and hindbrain), neural precursor proliferation, neural cell fate and migration, axon guidance, dendrite development and synapse formation.1 Given this broad spectrum of roles, it is no surprise that dysregulation of Wnt signaling causes a similarly broad spectrum of deleterious effects on neural development, and largely on this basis it has long been considered a candidate pathway in psychiatric pathogenesis. Its candidacy in this regard is supported by evidence that Wnt signaling pathway components are pharmacological targets of lithium,2, 3 hallucinogenic drugs 4 and antipsychotic medications,5, 6 as well as by evidence derived from behavioral tests in animal models7, 8, 9, 10, 11 and from positive associations in scattered human genetic studies.12, 13, 14, 15, 16, 17

Here we used massively parallel sequencing technology to simultaneously search for rare variants in 10 Wnt signaling pathway genes in a sample of Caucasian patients with ASD and healthy controls. We selected these 10 genes out of a much larger set of pathway molecules based primarily on empirical support for their individual relevance to psychiatry in human genetic (DKK1, DKK4, NEUROG1, CTNNB1),18, 19, 20, 21 animal behavior (DVL1, DIXDC1),7, 10 biochemical (WNT1, NEUROD1)22, 23 or other types of studies (PRICKLE4, DACT2).24, 25 Compared with whole exome sequencing strategies that are now underway, this very targeted approach benefited from our focus on functionally linked loci for which we have substantial expertise, allowing us to immediately follow up the most promising variant in in vitro functional assays.

Materials and methods

ASD samples (Caucasians)

Samples were lymphoblastoid cell lines from the Simons Simplex Collection. The diagnosis in these samples was based on the ADI-R (Autism Diagnostic Interview-Revised).26

Unscreened controls (Caucasians)

Healthy control subjects were recruited by the survey research company (Knowledge Networks, Menlo Park, CA, USA) from a nationally representative internet-based panel that was selected by random digit dialing.27 Subjects completed an online version of the CIDI-SF (Composite International Diagnostic Interview-Short Form) for lifetime history of anxiety, mood and substance use disorders. Subjects consented to anonymization and deposition of their clinical information and DNA in the National Institute of Mental Health (NIMH) repository for use in any medical research.

Pooling samples

Concentration of individual DNA samples was measured using Pico Green. Based on previous reports using an identical pooling strategy,28, 29 we expected sequencing error rates to be ∼0.5%. To distinguish a single real mutant allele in 876 chromosomes (minor allele frequency (MAF)=0.11%) from random sequencing artifacts, we pooled samples into 11 sets of 38–40 individuals each. Therefore, the frequency of a single mutant allele in each pool was 1%, significantly above the detection threshold of 0.5%.

Target enrichment

Each pool of genomic DNA from 38–40 individuals was enriched by PCR amplification for 10 genes involved in Wnt signaling (WNT1, DVL1, DIXDC1, PRICKLE4, DACT2, NEUROG1, DKK1, DKK4, NEUROD1, and CTNNB1), plus 1 control gene, DACT1, during the validation step (Supplementary Table S1). For each gene, the entire locus was amplified in 1 or 2 ‘long-range’ PCR reactions. The total amplification length for all 11 genes was 71.5 kb.

Sample preparation and high-throughput sequencing procedure

Each pool of target-enriched genomic DNA was used to build a prepped DNA library following Illumina sample preparation protocols. In brief, each pool was ligated with adapters containing a unique index allowing us to combine 11 pools in each sequencing lane of the HiSeq flow cell (HiSeq 2000, Illumina, San Diego, CA, USA). Finished flow cells carrying 438 individuals (1 lane containing 11 indexed pools of 38–40 individuals) were sequenced in the HiSeq 2000 for paired-end analysis.

Initial data processing (sequence assembly, SNP calling, quality assessment)

Bowtie30 was used for read alignment against the latest human genome reference. SAMtools,31 Picard (not published; http://picard.sourceforge.net/), BEDTools32 and the Genome Analysis Toolkit (GATK;33) were used for post-alignment processing. Multisample realignment around potential insertion/deletions (indels) and base quality score recalibration were performed before variant calling. For single-nucleotide polymorphism (SNP) calling, we used SyZyGy, an algorithm specifically designed to call bases in pooled samples.34

Analysis of conservation and prediction as deleterious

For each SNP, analysis by SeattleSeq (http://snp.gs.washington.edu/SeattleSeqAnnotation137/) provided information about base position conservation in the genome, whether it was intronic or exonic, and if exonic whether it was synonymous, missense or nonsense. Predicted deleteriousness for each SNP in exonic sequence was determined by PolyPhen-235 version v2.1.0r367 or more recent. If either version predicted that the SNP was deleterious, it was counted as deleterious.

Human cDNAs and mutagenesis

Human WNT1 (clone ID: 30915309) complementary DNA (cDNA) was purchased from Thermo Scientific (Rockford, IL, USA) and subcloned into a pcDNA 3.1(−) vector. The WNT1 S88R variant was generated using the QuikChange site-directed mutagenesis kit (Stratagene, La Jolla, CA, USA) and confirmed by Sanger sequencing. In subsequent functional assays, in order to avoid potential confounds arising from nonsystemic biases in measurements of plasmid yield or purity, each construct was prepared independently >2 times each by different experimenters in two different laboratories. Similar functional results were obtained from each independent preparation compared with wild type (WT).

Luciferase reporter assays

This was performed essentially as described previously.36 Briefly, both HEK293T and SH-SY5Y BAR-Luc cells were transfected with 1, 5, 20 and 50 ng pcDNA3.1(−) vector containing either WT WNT1 or WNT1-S88R. HEK293T cells were additionally transfected with 20 ng per well pBAR-Luc plasmid and 10 ng per well pRL-TK (Renilla luciferase control plasmid). Cells were cultured for 24 h, and then luciferase activity measured using the Dual-Luciferase Reporter Assay System (Promega, Madison, WI, USA) and an Envision multi-label plate reader (PerkinElmer, Waltham, MA, USA).

Quantitative reverse transcription-PCR

Total RNA from cultured HEK293T cells transfected with an empty pcDNA3.1(−) vector, or a pcDNA3.1(−) vector containing either WT WNT1 or WNT1-S88R, was isolated with the RNeasy Mini kit according to the manufacturer’s instructions (Qiagen, Hilden, Germany). The cDNA was synthesized from 2 μg of total RNA with Superscript III Reverse Transcriptase (Invitrogen, Carlsbad, CA, USA) using random hexamer primers. Quantitative reverse transcription PCR was performed using the PerfeCTa SYBR Green FastMix, ROX (Quanta, Byfleet, UK) with gene-specific primer pairs (WNT1 forward (Fwd): 5′-GCAGCGACAACATTGACTTC-3′, WNT1 reverse (Rev): 5′-GTGGCACTTGCACTCCTG-3′; AXIN2 Fwd: 5′-TTATGCTTTGCACTACGTCCCTCCA-3′; AXIN2 Rev: 5′-CGCAACATGGTCAACCCTCAGAC-3′) on an ABI 7900HT real-time PCR system (ABI Advanced Technologies, Columbia, MD, USA). Data were analyzed using the comparative CT method (User Bulletin No. 2, PerkinElmer Life Sciences, Boston, MA, USA). To normalize for loaded cDNA, glyceraldehyde 3-phosphate dehydrogenase (GAPDH) was used as an endogenous control (GAPDH Fwd: 5′-CAAGGTCATCCATGACAACTTTG-3′, GAPDH Rev: 5′-GGCCATCCACAGTCTTCTGG-3′). Results shown reflect combined data from two independent experiments, each done in triplicate reactions.

Statistical analyses

Unless stated otherwise in Figure legends, the data are expressed as mean±s.e.m. Data were analyzed by Student’s t-test, Fisher’s exact test or Pearson’s correlation test where appropriate, using GraphPad Prism software (GraphPad Software, La Jolla, CA, USA). PLINK software37 was used for Multidimensional Scaling and logistic regression. A P-value of <0.05 was considered to be significant.

Results

Using a next-generation (massively parallel Illumina-based) strategy, we sequenced promoter, 3' untranslated, intronic and coding regions of 10 Wnt pathway genes (WNT1, DVL1, DIXDC1, PRICKLE4, DACT2, NEUROG1, DKK1, DKK4, NEUROD1 and CTNNB1) in genomic DNA obtained from a first set of 198 Caucasian ASD patients and 240 healthy Caucasian controls, arranged in pools of 38–40 samples each. The number of reads was homogeneously distributed over the different amplicons and pools (and therefore between ASD and controls), reads mapped 88% of the targeted regions and coverage averaged 865-fold/individual. To validate the next-generation sequencing strategy employed, we compared results obtained using this strategy to results obtained using Sanger sequencing of the exons in a control locus (DACT1) in a subset of samples (40 ASD). Every DACT1 variant (n=8), including several singleton variants (n=3), found by Sanger sequencing in this set was also detected by the next-generation sequencing strategy, and the allele frequencies of all SNPs were identical using each methodology (Pearson’s correlation test; r=0.99; P<0.0001; Supplementary Table S1).

In the 10 genes under study, we identified 652 SNPs overall in this set of 438 combined ASD and control samples. We restricted most of our attention to nonsynonymous sequence variants found only in our own sample set (not referenced in dbSNP, 1000 Genomes, the Exome Variant Server (EVS) or the Exome Chip Design Consortium). To determine whether any of the putative SNPs found by this method were artifacts produced by the next-generation sequencing technology, we confirmed each of them by individual genotyping using either IPlex Sequenom or Sanger sequencing. Using this approach we confirmed 12 novel singleton missense variants distributed across 8 of the 10 loci included in our study (Supplementarty Table S2). Moreover, genotyping in a second set of 91 ASD and 144 control samples failed to identify any additional occurrences of each of these rare deleterious missense SNPs. Although each of these SNPs were found only once among the probands in our samples and had not been reported in any public database, genotyping of parents demonstrated that none of them had occurred de novo; all were transmitted from either the mother or father. Moreover, sequencing of family trios showed that each was also present in an unaffected sibling. This verifies that each of these variants is not an idiosyncratic artifact arising during lymphoblastoid cell line derivation, and suggests that they are extremely rare alleles (single-nucleotide variants) in the Caucasian population represented by our samples.

Although we found a greater number of extremely rare missense variants in our ASD sample compared with our control sample (4% in ASD vs 1.7% in controls; Figure 1), this difference in the ‘rare missense variant burden’ did not reach statistical significance (Fisher’s exact test, P=0.11). Nonetheless, when considering only those rare missense variants that were also predicted to be deleterious by PolyPhen-2, the difference in ASD versus control groups did become significant (3.5% in ASD vs 0.8% in controls; Fisher’s exact test, OR=4.37, P=0.04; Figure 1). We also examined the distribution of less informative variants (that is, synonymous, untranslated regions and intronic variants) in each group, and found a very similar distribution for all of them (synonymous: 3.5% in ASD vs 2.9% in controls; untranslated regions: 11.1% in ASD vs 9.6% in controls; intronic: 24% in ASD vs 21% in controls; Figure 1). The similarity in the incidence of less informative variants argues against the possibility that the difference in the incidence of rare deleterious variants can be attributed to differences in recent ancestry between the ASD and control sample groups.

Figure 1
figure 1

Burden of novel unique variants in the Wnt pathway in autism spectrum disorders (ASD). Percentage of individuals carrying a novel unique variant that was nonsynonymous and predicted to be deleterious by the PolyPhen-2 algorithm, compared with nonsynonymous, synonymous, untranslated region (UTR) or intronic single-nucleotide polymorphisms (SNPs). *P<0.05; NS, not significant.

Full size image

In addition to these novel SNPs, we found one previously referenced rare (<1%) missense SNP that was overrepresented in our ASD group versus our control group. This SNP (rs61758378) occurred in WNT1, the gene encoding the first-discovered mammalian Wnt ligand. To eliminate the possibility of a sequencing artifact we confirmed this result by genotyping all samples using IPlex Sequenom. We also genotyped an additional 91 ASD and 144 control samples to replicate and potentially extend these findings. In this way, we not only confirmed MAF results for rs61758378 obtained using the Illumina platform in the original sample set, we further found that the rs61758378 minor allele remained significantly overrepresented in the combined ASD sample compared with the combined control sample (8 A/T in 267 ASD (MAF=1.69%) vs 1 A/T in 377 controls (MAF=0.13%); Fisher’s exact test, allelic P=0.0048; Table 1). Furthermore, when we compared the rs61758378 MAF in our combined ASD sample with the reported MAF in the EVS for a much larger group of individuals, the association also remained significant (8 A/T in 267 in our combined ASD sample (MAF=1.69%) vs 45 A/T in 4300 European/American individuals in the EVS (MAF=0.5%), Fisher’s exact test, allelic P=0.011). Because different allele frequencies are reported in European populations (TSI, Tuscans in Italy, MAF=3.0%; CEU, Northern Europeans in Utah, MAF=0.5%), we wanted to evaluate and address possible effects of subtle stratification in ancestry between cases and controls. We therefore performed Multidimensional Scaling using genome-wide SNP genotype data for a subset of the total samples used in this analysis (235 SSCs and 329 controls) using PLINK. Multidimensional Scaling using reference CEU and TSI data showed that the ratio of individuals clustering with CEU and TSI was not significantly different between ASD and controls. Nevertheless, we noted that one of the first four dimensions was significantly different between ASD and controls, indicating possible stratification. We sought to account for this using a logistic regression analysis including the first four dimensions as covariates. In this analysis, the odds ratio (OR) remained high (OR=8.17) and the P-value remained close to significance (P=0.053). We note that our analyses have low power for SNPs with such a low MAF and in a modest number of samples.

Table 1 Association of WNT1 rs61758378 with ASD
Full size table

The rs61758378 minor allele changes a conserved serine at position 88 into arginine. Although this change is not predicted to be deleterious by PolyPhen-2 (HumDiv score=0.174), serine 88 is nonetheless a very highly conserved residue across species (PhasCons score=0.998), and the change to arginine at this residue is predicted to be deleterious by Grantham score, which categorizes codon replacements into classes of increasing chemical dissimilarity (Grantham score=110, moderately radical change).38 Given the potential association between the WNT1 rs61758378 minor allele and ASD in our sample set, we therefore tested for functional effects of this variant on Wnt signaling pathway activation. To accomplish this, we engineered this variant into a human WNT1 cDNA (WNT1-S88R) using site-directed mutagenesis and then tested it alongside WT Wnt1 in several assays for Wnt signaling activity.

First, we performed a plasmid-based Wnt reporter assay in human embryonic kidney (HEK293T) tissue culture cells. This is a well-established method to quantify transcriptional activation downstream of the Wnt/β-catenin signaling pathway, based in this case on transient transfection of the pBAR-Luc reporter plasmid that expresses luciferase downstream of a β-catenin-responsive promoter.36 Cells were co-transfected with pBAR-Luc and a dose range (1–50 ng) of plasmids encoding WNT1-S88R or WT WNT1 to test for the ability of these proteins to activate Wnt/β-catenin pathway-dependent transcription. At all doses tested, WNT1-S88R showed increased activity in this assay compared with WT WNT1 (Figure 2a). Moreover, this difference was statistically significant at the 20 ng dose (WT WNT1, 1.00±0.09 vs WNT1-S88R, 1.37±0.087, t-test, P=0.02) and at the 50 ng dose (WT WNT1, 1.00±0.18 vs WNT1-S88R, 1.86±0.22, t-test, P=0.016; Figure 2a). Similar results were obtained when experiments were replicated in a different laboratory by a different investigator using independently prepared DNA constructs and reagents (pSuperTop; data not shown).

Figure 2
figure 2

Functional analysis of S88R. (a) Plasmid-based Wnt/β-catenin pathway reporter assay in HEK293T cells. HEK293T cells were transfected with the pBAR-Luc reporter plasmid together with plasmids encoding either WNT1-S88R or wild-type (WT) WNT1 at doses ranging from 1 to 50 ng. (b) Quantitative reverse transcription-PCR (RT-PCR) of the Wnt/β-catenin pathway endogenous target gene AXIN2 in HEK293T cells transfected with 5 ng of each WNT1 plasmid, expressed as a ratio to WNT1 expression also determined by RT-PCR. (c) Luciferase reporter activation in SH-SY5Y BAR-Luc cells transfected with 50 ng of each WNT1 plasmid. *P<0.05 and ***P<0.001.

Full size image

The aforementioned assays rely on transient transfection of reporter plasmid, and on this basis they are conceptually more variable than assays based on read-out of an endogenous Wnt pathway target gene or of a single-copy Wnt pathway reporter integrated into the genome. Accordingly, we confirmed these results in the same (HEK293T) cell line using quantitative reverse transcription-PCR to directly measure transcriptional activation of the endogenous Wnt/β-catenin pathway target gene AXIN2.39 Mirroring results obtained with the Wnt reporter assay, compared with WT WNT1, WNT1-S88R showed significantly increased transcriptional activation of AXIN2 (WT Wnt1, 1.00±0.025 vs Wnt1-S88R, 2.16±0.43, t-test, P=0.02; Figure 2b). We confirmed these results via a third in vitro assay that relies on SH-SY5Y BAR-Luc cells. This is an independently derived immortalized cell line with a single copy of a Wnt/β-catenin pathway luciferase reporter stably integrated into its genome.36 Consistent with data obtained by the plasmid-based Wnt reporter assays and by quantitative reverse transcription-PCR in HEK293T cells, SH-SY5Y BAR-Luc cells showed significantly increased activation of their reporter after transfection of WNT1-S88R (WT Wnt1, 1.00±0.04 vs Wnt1-S88R, 1.28±0.03, t-test, P=0.0003; Figure 2c).

Taken together, our data from cell-based transcription target assays suggest that the rs61758378 minor allele associated with ASD in our sample (WNT1-S88R) encodes a modestly more active form of the WNT1 ligand compared with the most prevalent WNT1 allele found in Caucasians.

Discussion

Our finding in ASD of a higher burden of rare missense variants distributed across 7 of 10 Wnt signaling pathway genes tested, and of an ASD-associated functional variant at the WNT1 locus, supports that dysfunction of this pathway contributes to ASD susceptibility.

In functional assays for the Wnt/β-catenin pathway, the biochemical cascade physiologically activated by the WNT1 gene product,40, 41 WNT1-S88R consistently showed increased activity. In animal models, hyperactivation of Wnt/β-catenin signaling can lead to dominant neurodevelopmental defects including in brain regionalization,42, 43 neural precursor proliferation, cell fate specification, neuron migration,44, 45, 46, 47 and in dendrite and synapse formation or function.48, 49 A mutation in a gene, such as WNT1, that modestly increases activation of this downstream signaling pathway would therefore be expected to result in neurodevelopmental and behavioral phenotypes. Moreover, as ASD and other major psychiatric disorders are generally not associated with gross neuroanatomical defects, changes in Wnt signaling activity contributing to susceptibility for these disorders should be modest. Anatomical changes resulting from a signaling variant such as WNT1-S88R might only be evident at a cellular level (for example, in neural subtype numbers and position, neural morphology, connectivity and so on) or at a subcellular level (for example, in dendritic spine morphology, synapse numbers, synapse strength, synapse plasticity and so on), yet could still be sufficient to cause behavioral symptoms. This of course conforms closely to the central features of major psychiatric disorders, which are distinguished from neurological disorders in part by the absence of focal or grossly evident structural brain defects. Moreover, it is not our intent to posit that this or the other rare coding variants we have identified here in WNT1 and in Wnt signaling pathway component genes are monogenic causes of psychiatric illness.50 Genetic susceptibility to these disorders is highly heterogeneous and thought to be frequently multigenic; that is, for most patients genetic loading is thought to be distributed across many loci that affect a common pathogenic pathway.51, 52 In this susceptibility model, an allelic variant at any one locus may only appreciably increase disease risk in the context of variants at other loci that affect the same pathway in the same individual. The rare nonsynonymous sequence variants found in both affected and unaffected individuals within the same family in this study may contribute to ASD in the context of other genetic and nongenetic factors that increase susceptibility.

Several lines of evidence support that Wnt ligands (of which there are 19 family members in humans) participate in ASD susceptibility. A recent whole exome sequencing study identified copy number variants affecting both WNT3 and WNT9B in an ASD patient,53 whereas another study showed that WNT3 expression is increased in the prefrontal cortex of young autistic patients compared with controls.54 Rare variants in WNT2 have been reported to segregate with ASD.16 Although other groups have not always replicated this linkage,55, 56 two more recent studies have found SNPs in WNT2 associated with either ASD or speech delay.15, 17 In a Valproate-exposure rat model of ASD, expression of both WNT1 and WNT2 was upregulated in the developing forebrain and this was associated with increased Wnt/β-catenin pathway activity in that brain region.57 Other recent sequencing studies have identified several ASD-susceptibility genes that encode functionally linked components of a protein network involved in the regulation of β-catenin-dependent transcription, as occurs in Wnt/β-catenin signaling.21, 58, 59 The findings of the present study add to this evidence by suggesting that minor alleles and extremely rare variants at the WNT1 locus and other Wnt pathway loci may contribute to risk for this psychiatric disorder, a possibility that merits further attention in ongoing human genetic studies focused on ASD. Given evidence that genetic etiologies are shared across several psychiatric disorders,60, 61, 62 as well as that Wnt signaling is involved more broadly in psychiatric pathogenic and therapeutic mechanisms,63 these loci should also be carefully scrutinized in ongoing human genetic studies of other major psychiatric disorders, such as schizophrenia and major affective disorders.