Novel Autism Subtype-Dependent Genetic Variants Are Revealed by Quantitative Trait and Subphenotype Association Analyses of Published GWAS Data

Quantitative Trait Association Analyses

The flowchart in Fig. 1 describes the experimental design and analyses that were used to derive the combined set of 18 novel and statistically significant SNPs that associate with 4 separate subtypes of ASDs. Raw item scores from the ADI-R score sheets of 2939 ASD cases were summed for spoken language skills, non-verbal communication, play skills, social development, and insistence on sameness/rituals, as described previously [13]. The specific items used to obtain the total score per “trait” category for each individual are shown in Table S1 and the profiles of total scores for each category are shown for the 2939 individuals in Fig. S1. We then conducted quantitative trait association analyses using the distribution of scores in each of the categories to identify sets of SNPs that associate with symptomatic severity of each of the behaviors listed in Table S1. A nominal p-value≤10⁻⁵ was used in this “discovery phase” designed to screen candidate SNPs that are potentially of functional significance with respect to the ASD traits. These sets of symptom-associated SNPs (or quantitative trait loci, QTL) are shown in Table S2.

ASD subtype-dependent genetic association analyses with trait-associated SNPs

Next, we performed cluster analyses as described by Hu and Steinberg [13] to divide the autistic cases into 4 non-overlapping phenotypic subgroups according to symptomatic severity profiles derived from 123 items on the ADI-R assessments. We postulated that this subtyping procedure, which reduces the behavioral/symptomatic heterogeneity among the cases within each subgroup, is likely to also restrict the genetic heterogeneity within each subgroup. The results of K-means cluster analyses (K = 4) of the ADIR data from the 2939 individuals (which included 1867 genotyped cases from the GWA study) is shown in Fig. S2A. An unsupervised principal components analysis (PCA) was used to demonstrate that the 4 subgroups of individuals with ASDs from the K-means cluster analyses results in effective separation of cases according to similarity of symptomatic profile as defined by the 123 item scores from the ADI-R (Fig. S2B). The resulting phenotypic subgroups were then used in genetic association analyses with each set of QTL derived from the five quantitative trait association analyses, where the cases were either divided into the 4 ASD subtypes or used as a combined case group against the 2438 nonautistic controls. These analyses produced 5 sets of QTL that were associated with specific subtypes of ASDs ( Tables 1 , 2 , 3 , 4 , 5 ). It is interesting to note that SNPs associated with certain genes (such as HTR4, BACH2, and TINAG among Language QTL, and OR5B3, OR5B17, and ERBB4 among Nonverbal QTL) are replicated in more than one ASD subtype. Finally, significant SNPs with Bonferroni-adjusted p-values≤0.05 from each of the 4 separate subtype-dependent association analyses with QTL were combined into a single set containing 18 unique SNPs, and case-control association analyses were performed for which the cases were either combined into one group or divided according to subphenotype.

Partial replication of SNPs between subtypes of ASDs

Table 6 shows the SNPs associated with each subtype of ASDs that resulted from the final association analyses using the combined highly significant QTL and 4 non-overlapping subgroups of ASD cases. Eighteen of the SNPs have p-values<0.05 in at least one subtype even after using the stringent Bonferroni correction for multiple comparisons. Note that 10 of the SNPs, including rs317985, rs7785107, rs11671930, rs7950390, rs12266938, rs3861787, rs7725785, rs1827924, rs1231339, and rs757099, are associated with more than one subtype. Two of the replicated SNPs (rs317985 and rs7785107) are significant in two subtypes after Bonferroni adjustment (p<0.05), while the remaining 8 exhibit lower levels of significance (nominal p-values from 0.0037–0.051 or FDR_BH adjusted p-values of 0.0087–0.088) in the second (or third) subtype. Association of these QTL with more than one subtype of ASD serves as a replication for these 10 SNPs. While the replication of SNPs in different subtypes may be interpreted as evidence for those SNPs being highly penetrant for ASDs, the subtype-dependent differences in minor allele frequency (MAF) and odds ratios (OR) associated with the shared SNPs suggest that the subtypes are genetically heterogeneous. Figure 2 summarizes the extent of SNP overlap among the 4 ASD subtypes (case cohorts) and clearly demonstrates that the odds ratios are distinctly different for different subtypes that share the same SNP. All of the QTL associated with specific genes are present in non-exonic (promoter, intronic, or intergenic) regions. Interestingly, all but one of the SNPs residing within intergenic regions can be associated by band position to rare copy number variants (CNV) that have been recently identified for ASDs [15], thus providing further support for the probable relevance of these SNPs to ASDs. These CNVs are noted in Table 6.

Effect of subtyping of ASDs on association of previously identified SNPs within Chr5p14.1

Because none of the SNPs identified in this re-analysis overlapped with those of the previously published genome-wide association study [9], we examined the association of the 6 SNPs that were reported in the published study with our subphenotypes. Table 7 shows that only the “Moderate” ASD subtype (363 cases) is associated with two of the SNPs, with Bonferroni-adjusted p-values of 0.035 and 0.053. Interestingly, these 2 SNPs have the least significant combined p-values in the published study. The remaining 4 SNPs were suggestively significant with FDR_BH-adjusted p-values of 0.074 in this subtype. The combined cases (1867 individuals in all) as well as the other 3 ASD subtypes show no association with any of the 6 SNPs even though there are more cases in each of these groups than in the Moderate group. This finding further illustrates the value of analyzing subphenotypes of ASDs in genome-wide association analyses in order to reduce the genetic heterogeneity of the population studied.

Pathway analyses of SNP-containing genes

To obtain a better understanding of how the novel SNPs identified in this study potentially relate to the biology of autism, pathway analysis was conducted to develop a better sense of the relationships among the SNP-associated genes and their impact on higher level functions and diseases. Fig. 3 shows a gene network constructed using Pathway Studio 7 which includes seven of the 9 genes associated with SNPs found within gene promoters or introns. Of the 7 genes, HTR4 and GCH1 are “hubs” showing the highest “connectivity” with other components within the network. The specific relationships between these two genes and other network components are more clearly illustrated in Figures S3 and S4. It is noted that many of the cellular and higher level processes in this network, such as neurogenesis, axonogenesis, steroid metabolism, cell proliferation, long-term synaptic potentiation, learning and memory are relevant to identified deficits in ASDs.

Quantitative trait and phenotype incorporation into genetic analyses

There have been relatively few large-scale genetic studies on ASDs that incorporate quantitative traits or phenotypes into genome-wide analyses. The value of using phenotypic subsets of ASDs in genetic analyses was demonstrated in an early linkage study by Nurmi et al. who found that subsetting the subjects (and families) according to presence or absence of savant skills increased the HLOD score of the D15S511 marker within the GABRB3 gene on chr15q11–q13 from 0.8 in the entire study population to 2.6 for the savant-skills-positive subset of individuals [16]. Later, Alarcon et al. used quantitative trait analysis of 2 language scores on the ADI-R combined with ordered subsets analysis for subject stratification to reveal a strong language QTL on 7q35, and observed additional linkage peaks on several other chromosomes without stratification [17]. Chen et al. used a combination of 7 ADI-R scores describing nonverbal communication and the same approach to uncover QTL for nonverbal communication [18]. Similarly, Duvall et al. used quantitative trait analysis based upon scores from the Social Responsiveness Scale to discover QTL for social endophenotypes with strong linkage to chromosomes 11 and 17 across all samples and additional linkages to chromosomes 4, 8, and 10 in “male-only” families [19]. More recently, Nijmeijer et al. identified loci overlapping ASDs and attention-deficit/hyperactivity disorder (ADHD) in a genome-wide study of children with ADHD and autistic traits (but not classical or atypical autism) by using total and subscale scores from the Social Communication Questionnaire (an autism diagnostic instrument) as well as by including scores related to severity of ADHD as covariates [20]. Along a similar vein, Ronald et al. examined both social and non-social autistic-like traits in the general population, quantified on the basis of DSM-IV social and non-social scales, to tease out 22 and 20 SNPs that associated with these traits, respectively [21]. However, only one SNP survived testing in an independent sample and its p-value was still only nominally significant in an all male sample. Recently, quantitative trait analysis was applied in a genome-wide linkage study of two repetitive ASD behavior phenotypes, “insistence on sameness” (IS) and “repetitive and stereotyped behavior” (RSMA) [22]. In this study, Cannon et al. demonstrated IS and RSMA were linked to 2 separate loci, each exhibiting respective HLOD scores of 3.42 and 3.93, suggesting that more narrowly-defined phenotypes improve linkage signals. When social behavior, as reflected by severity scores on the Social Responsiveness Scale was used as a quantitative trait, Coon et al. not only replicated a susceptibility locus identified previously in the Duvall et al. study [19], but also identified additional linkage peaks with HLOD scores ranging from 2.23–3.64 [23]. Interestingly, St. Pourcain et al. recently reported that one of the SNPs identified at Chr5p14 by GWA analyses [9] also associates by quantitative trait analyses with the social communication phenotype [24]. They suggested that this SNP, rs4307059, may also associate with multiple subthreshold impairments across a variety of ASD-related functions.

In the present study, we observe that rs4307059 is suggestively associated only with the Moderate ASD subtype, with an unadjusted p-value of 0.074 (Table 7). Thus, the above-mentioned studies, in combination with ours, reinforce the concept that the incorporation of quantitative traits into genetic analyses greatly improves the ability to elicit significant genetic information on complex disorders such as ASDs. Furthermore, the genes associated with the QTL in Tables 1, 2, 3,4, 5 and Table S2 may provide clues to the biology underlying the various autistic traits. Our study, which uniquely combines QTL with ASD subtype-dependent genetic association analyses, provides an additional refinement to genome-wide studies that results not only in the identification of a greater number of significant SNPs, but also concomitant and internal validation of SNPs that are replicated between and among the ASD subtypes which represent independent cohorts of affected individuals. Since the number of individuals in each subgroup ranges from 363–639, our study also demonstrates the increased statistical power inherent in this approach. Clearly, there are many ways that individuals with ASDs can be divided into subgroups (such as by regression following normal development, co-morbid symptoms, IQ, etc.), and some of these methods may prove to be superior to the method described here and in our previous study [13]. However, the importance of using phenotypically distinct cohorts in large-scale genetic association analyses is illustrated here using a method that is capable of distinguishing meaningful phenotypes across the whole spectrum of ASDs (severe to mild), based on severity profiles across a broad range of behaviors and symptoms characteristic of ASDs.

Biological relevance of SNP-containing genes

To examine the biological processes and pathways that might be impacted by the SNP-associated genes, we performed pathway analyses on the genes using Pathway Studio 7 software. Fig. 3 shows that 7 of the 9 SNP-containing genes are included in a gene network with many genes leading to “axonogenesis,” a function which has been reported to be enriched among genes overlapping copy number variants in autistic patients [25]. HTR4 and GCH1 are network “hubs” connecting with many other genes, cellular processes and disorders (see Figs. S3 and S4 for specific connections). As illustrated in Fig. S3, HTR4 [5-hydroxytryptamine (serotonin) receptor 4] regulates neurogenesis, long-term synaptic potentiation and, in turn, learning and memory, as well as the release of neurotransmitters (dopamine, acetylcholine), peptide hormones (AVP, OXT, PRL, VIP) and steroid compounds (cortisol, corticosterone). Thus, any alteration in the expression or function of this gene can be expected to have wide-ranging consequences on processes known to be affected by ASDs. It is notable that one of the SNPs associated with HTR4, rs7725785, is associated with three ASD subtypes (Table 6). However, the odds ratio is 1.44 for the severe language-impaired subtype while it is 0.68 and 0.74 for the moderate and mild subtypes, respectively. Interestingly, a reduction in HTR4 expression was observed only for the language-impaired subtype of ASDs and not for the moderate (savant) or mild subtypes in our earlier gene expression study [14], suggesting that this SNP, which is located at an intron boundary within HTR4, may play a role in level of expression. Genetic variants in HTR4 have also been associated with schizophrenia [26], bipolar disorder [27] and attention deficit/hyperactivity disorder [28], further suggesting some functional overlap among these neurological disorders. More recently, a de novo translocation on chromosome 5 close to HTR4 has been identified in an autistic boy [29], lending additional support for the relevance of this gene in ASDs.

The other hub gene, GCH1 [GTP cyclohydrolase I], is the rate-limiting enzyme in the de novo biosynthesis of tetrahydrobiopterin which is, in turn, required for the biosynthesis of folate, serotonin, dopamine, and catecholamines (Fig. S4), all of which are important for neural development and functions. It is interesting to note that elevated expression of GCH1 has been implicated in mood disorders [30], while genetic polymorphisms or mutations in GCH1 have been associated with pain sensitivity [31]–[33], and dystonia [34], which are often associated with ASDs. Although these genes are not likely to be causal for ASDs, genetic polymorphisms in them may be associated with some of the comorbid symptoms or pathobiology of ASDs.

Together, these analyses provide support for the possible biological relevance of these QTL to ASDs and identify additional candidate genes, such as HTR4 and GCH1, for functional testing. Importantly, this study also reveals genetic biomarkers which not only may be used for diagnostic screening of ASDs, but also offer the additional advantage of being associated with ASD subtypes that may be linked to specific and targeted therapies through pharmacogenomics studies. Finally, the association of different SNPs with the 4 subtypes of ASDs described here reinforces the idea that there are multiple genetic etiologies giving rise to the autistic spectrum, while the shared SNPs between different subtypes may reveal common genetic mechanisms responsible for core symptoms.

Summary

This study is the first to demonstrate the value of using a combination of quantitative trait analysis and subphenotyping of individuals with ASDs to identify genetic variants (SNPs) that associate with specific behavioral phenotypes of ASDs. It is noted that no Bonferroni-significant SNPs are detected when all 1867 autistic cases are combined into a single group and compared against 2438 non-autistic controls. Thus, even though the number of cases is lower in each of the subgroups, there is more power to detect statistically significant SNPs associated with the more homogeneous subgroups of individuals with ASDs than with the combined group of cases. Subtyping also creates independent case cohorts which are shown to replicate 10 of the novel SNPs identified in this study, thus providing a form of internal validation. Differences in minor allele frequencies of these 10 SNPs in the different cohorts further demonstrate the genetic heterogeneity among the subtypes. Together, these findings not only reveal novel subtype-dependent candidate genes for ASDs, but also identify potential genetic markers for diagnostic screening.

GWA data used for this study

Genome-wide association data from the study by Wang et al. [9] were downloaded from the Autism Speaks website at ftp://ftp.autismspeaks.org/Data/CHOP_PLINK/AGRERELEASE.ped. For this study, we used the file named CHOP.clean100121 where the data were “cleaned” by Jennifer K. Lowe in the laboratory of Daniel H. Geschwind, M.D., Ph.D. at UCLA. The cleaning procedure involved extensive sample and pedigree validation, exclusion of SNPs a) missing >10% data, b) with HWE p<0.001, c) with MAF<0.01, and d) with >10 Mendelian errors. The final dataset included 4327 genotyped individuals and 513,312 SNPs on the Illumina HumanHap550 Bead Chip.

ADI-R data

ADI-R assessments for 2939 individuals were obtained from Autism Speaks through Dr. Vlad Kustanovich of the Autism Genetics Resource Exchange. Of these, 1867 individuals were among the cases genotyped by Wang et al. [9]. Scores of 123 items on the ADI-R score sheets of each individual were analyzed as described by Hu and Steinberg [13] to identify ASD subtypes (that is, phenotypic subgroups) which are represented in Fig. S2.

Determination of quantitative traits

Raw item scores from the Autism Diagnostic Interview-Revised (ADI-R) score sheets of 2939 ASD cases were summed for spoken language skills, non-verbal communication, play skills, social development, and insistence on sameness/ritualistic behaviors. The specific items used to obtain the total score per category for each individual, shown in Table S1, were noted in an earlier study [13] to exhibit average differences in severity among several subtypes of ASDs, described below. The sums of item scores within each of the 5 categories were used as quantitative traits for genetic association analyses using the genotype data reported by Wang et al. [9]. Profiles of the traits across the 2939 individuals are shown in Fig. S1.

Subtyping of individuals with ASDs

Two thousand nine hundred and thirty nine (2939) individuals with ASDs were divided into phenotypic subgroups using clustering tools within the MeV software package [35], as previously described by Hu and Steinberg [13]. Briefly, subtyping of individuals with ASDs involved K-means cluster analysis (with K = 4) of scores from 123 items from the ADI-R score sheets of each individual which were adjusted as described to fall into a severity range of 0 (normal) to 3 (highest severity). A figure of merit analysis (not shown) indicated that the individuals with ASDs were optimally represented by 4 subgroups. A principal components analysis (PCA) was then used to verify that the 4 subgroups of individuals identified by the K-means cluster analyses were distinguishable by this unsupervised test. Fig. S2 shows the symptomatic profiles of the 4 ASD subtypes as well as their separation into discernible clusters by PCA. The subtypes are named “Language-impaired,” “Intermediate,” “Moderate,” and “Mild,” and contain 639, 478, 363, and 387 cases, respectively.

Quantitative trait association analyses

Using the score totals for the 5 categories of autistic symptoms exhibited by each of 1867 cases as quantitative traits, we utilized PLINK [36], a program to analyze whole genome association data, to perform quantitative trait association analyses with the genotype data reported by Wang et al. [9]. From the results of each of the 5 analyses (Table S2) in this discovery phase, we selected SNPs with unadjusted p-values<10⁻⁵, which prioritized SNPs filtered by association with quantitative traits relevant to ASDs. These filtered sets of SNPs were then used in case-control association analyses, as described below.

Case-control association tests

Using the 5 sets of QTL generated for each symptom category, association analyses using PLINK were performed with each set where the cases, in addition to being treated as a single group of 1867 cases vs. 2438 controls, were also divided into the 4 ASD subtypes that were determined by the ADI-R cluster analyses described above. It should be noted that each ASD subgroup represents an entirely separate cohort of cases. From each of the 5 sets of genetic association analyses with subtypes and QTL ( Tables 1 , 2 , 3 , 4 , 5 ), we selected SNPs associated with each ASD subtype with a Bonferroni-adjusted p-value≤0.05 and combined them (a total of 18 unique SNPs) for a second case-control association analysis using the combined and subtyped ASD cases against controls. The results of this final analysis are shown in Table 6 .

Pathway analysis

Pathway Studio 7 software (Ariadne Genomics, Inc.) was used to generate relational gene networks using the SNP-containing genes listed in Table 6.