Association of oxytocin receptor (OXTR) gene variants with multiple phenotype domains of autism spectrum disorder

Autism spectrum disorder (ASD) is characterized by core deficits in social behavior, communication, and behavioral flexibility. Several lines of evidence indicate that oxytocin, signaling through its receptor (OXTR), is important in a wide range of social behaviors. In attempts to determine whether genetic variations in the oxytocin signaling system contribute to ASD susceptibility, seven recent reports indicated association of common genetic polymorphisms in the OXTR gene with ASD. Each involved relatively small sample sizes (57 to 436 families) and, where it was examined, failed to identify association of OXTR polymorphisms with measures of social behavior in individuals with ASD. We report genetic association analysis of 25 markers spanning the OXTR locus in 1,238 pedigrees including 2,333 individuals with ASD. Association of three markers previously implicated in ASD susceptibility, rs2268493 (P = 0.043), rs1042778 (P = 0.037), and rs7632287 (P = 0.016), was observed. Further, these genetic markers were associated with multiple core ASD phenotypes, including social domain dysfunction, measured by standardized instruments used to diagnose and describe ASD. The data suggest association of OXTR genetic polymorphisms with ASD, although the results should be interpreted with caution because none of the significant associations would survive appropriate correction for multiple comparisons. However, the current findings of association in a large independent cohort are consistent with previous results, and the biological plausibility of participation of the oxytocin signaling system in modulating social disruptions characteristic of ASD, suggest that functional polymorphisms of OXTR may contribute to ASD risk in a subset of families.


Introduction
Autism spectrum disorder (ASD) is characterized by abnormalities in three domains: social interaction deficits, language impairments, and repetitive behaviors with restricted interests. Despite a widely varying clinical presentation, ASD is highly heritable, with studies demonstrating a concordance rate of 80-95% in monozygotic twins and 0-31% in dizygotic twins (Folstein and Rutter 1977;Ritvo et al. 1985;Bailey et al. 1995;Le Couteur et al. 1996;Taniai et al. 2008;Rosenberg et al. 2009). However, genetic linkage studies have indicated signals for linkage on nearly every chromosome (International Molecular Genetic Study of Autism Consortium 1998; Philippe et al. 1999; Barrett et al. 1999;Risch et al. 1999;International Molecular Genetic Study of Autism Consortium 2001;Yonan et al. 2003;Ylisaukko-oja et al. 2006), suggesting that multiple genes contribute to ASD susceptibility.
The gene encoding the oxytocin receptor, OXTR, is a strong functional ASD candidate gene based on its known role in modulating social behavior (Ebstein et al. 2009;Hammock and Levitt 2006). Pharmacological and genetic manipulations have demonstrated a causal role for oxytocin (OXT) and its receptor (OXTR) in the regulation of species-typical social behavior. OXT facilitates social recognition behavior (Ferguson et al. 2000) and modulates maternal behavior (Pedersen and Prange 1979;Kendrick et al. 1987). Moreover, OXT signaling facilitates social preferences between adult monogamous rodents (Williams et al. 1992;Williams et al. 1994). In humans, there have been a number of recent studies demonstrating enhanced functions relevant to social behavior following oxytocin application in healthy adults (Heinrichs et al. 2003;Kosfeld et al. 2005;Rimmele et al. 2009;Zak et al. 2007;Domes et al. 2007a;Domes et al. 2007b;Guastella et al. 2008a;Guastella et al. 2008b;Hurlemann et al. 2010). Studies probing the relationship between OXT and ASD have generated complex findings. Reported reductions in plasma OXT (Green et al. 2001;Modahl et al. 1998) and increases in unprocessed OXT peptides (Green et al. 2001) in young children with ASD are contrasted by measures of higher baseline levels of serum OXT in young adults with ASD compared to control subjects (Jansen et al. 2006). Intravenous or intranasal applications of OXT in youth (Guastella et al. 2010) and adult (Hollander et al. 2007;Hollander et al. 2003;Andari et al. 2010) cohorts with ASD generally improved various functions relevant to social behavior.
Genetic analyses of OXTR contribution to ASD risk have been inconsistent. Most genome-wide linkage analyses did not find a peak near the OXTR gene on chromosome 3p25 (International Molecular Genetic Study of Autism Consortium 1998; Philippe et al. 1999; Barrett et al. 1999;Risch et al. 1999;International Molecular Genetic Study of Autism Consortium 2001). However, one of the largest linkage studies performed to date, including 314 families, highlighted a linkage peak directly over the OXTR gene, establishing OXTR as a positional candidate gene (Ylisaukko-oja et al. 2006). Seven genetic association studies of OXTR with ASD have been reported, but consistent association of alleles is lacking. A study of 195 Chinese families indicated association of two markers in the third intron of OXTR, rs2254298 A allele and rs53576 A allele (Wu et al. 2005). Three family-based studies failed to find association of rs2254298 (Liu et al. 2010;Lerer et al. 2008;Wermter et al. 2010) and another study did not find association of rs53576 (Jacob et al. 2007). However, evidence for association of the rs2254298 A allele was found in a Japanese case-control sample (Liu et al. 2010). In Caucasian families, the opposite allele, the rs2254298 G allele, was associated with ASD risk in one study (Jacob et al. 2007) and contributed to a risk haplotype in another study (Lerer et al. 2008). A third study of Caucasian families did not test rs2254298 directly, but found association of a nearby intron 3 marker, rs2268493 (Yrigollen et al. 2008). In addition to the evidence for association of OXTR intron 3 markers, markers near the 3′ untranslated region (UTR) have also been implicated in ASD risk (Lerer et al. 2008;Tansey et al. 2010). Lerer et al. (2008) described association of the 3′ UTR marker rs1042778 G allele and Tansey et al. (2010) found association of the intergenic marker rs7632287 with ASD risk. While different reports indicate positive association of OXTR markers in intron 3 and the 3′ UTR, the studies involved small family cohorts of less than 500 families each, increasing the risk for spurious associations (Sullivan 2007). Further, in those studies that included analysis of phenotype data, the OXTR markers associated with ASD failed to show association with measures of social behavior in individuals with ASD (Lerer et al. 2008;Yrigollen et al. 2008).
Given the sound biological and mixed genetic evidence in favor of OXTR contribution to ASD, we hypothesized that analysis of a large sample would provide greater power to detect association of OXTR markers with ASD susceptibility. Here, we report association analysis of 25 genetic markers spanning OXTR in 1,238 families, including 2,333 individuals with ASD. In addition to analysis of association with ASD diagnosis, we examined association with quantitative scores derived from three instruments used to diagnose and describe autism phenotypes: the Autism Diagnostic Interview-Revised (ADI-R), the Autism Diagnostic Observation Schedule (ADOS), and the Social Responsiveness Scale (SRS).

Methods and materials
Sample The family-based sample consisted of 5,432 individuals from 1,238 pedigrees and included 2,333 individuals with ASD (Table 1). The majority of the sample (921 pedigrees) was collected by the Autism Genetics Resource Exchange (AGRE) Consortium. The remaining 317 "non-AGRE" pedigrees were collected at the University of Iowa, Stanford University, Tufts University, and Vanderbilt University. Approximately 91% of the families had more than one child with ASD (multiplex). The genotyped sample is 95% Caucasian by self-report. The sample analyzed here does not overlap any of the previous reports of genetic association studies for OXTR in ASD (Wu et al. 2005;Liu et al. 2010;Lerer et al. 2008;Wermter et al. 2010;Jacob et al. 2007;Yrigollen et al. 2008;Tansey et al. 2010). Scores on the ADI-R, ADOS, and SRS were available for 2,000 individuals from AGRE families ( Table 2). All research was approved by the Institutional Review Boards of Vanderbilt University and the University of Southern California.
In-house SNP genotyping and quality control Genotyping was performed using TaqMan ™ SNP genotyping assays on the ABI Prism 7900HT and analyzed with SDS software as previously described (Campbell et al. 2006;Campbell et al. 2008). SNP Genotyping Assays-On-Demand were obtained from Applied Biosystems (Foster City, CA). Genotyping was performed in a 384-well plate format using 3 ng genomic DNA. Quality control measures included seeding of each 384-well plate with eight to ten blank negative control wells and 20-30 duplicated positive control samples. Automated allele calls were made with SDS Data Collection software and reviewed by an experienced operator according to protocol. The overall no-call rate was <5% for each of the assays. All analyzed markers were in Hardy-Weinberg equilibrium (HWE; P>0.05).
GWAS genotyping GWAS of the AGRE Consortium sample was conducted by the Broad Institute on the Affymetrix 5.0 platform, which includes over 500,000 SNPs, and made available publicly. We downloaded the genotype information from the AGRE website (www.agre.org). For analysis purposes, we used Whole-genome Association Study Pipeline. We set the marker genotyping efficiency threshold to a minimum of 95%, the minor allele frequency threshold to a minimum of 0.05, the Hardy-Weinberg equilibrium threshold to a minimum of 0.01, and the marker Mendelian error rate to a maximum threshold of 15 errors. Any marker that did not meet any of the preceding specifications was removed from all further analyses. In addition, all monomorphic markers and all markers composed of only heterozygous genotypes were removed. Following this whole-genome quality control, we selected the markers in the OXTR gene plus 10 kb of flanking region on each side (a total of 39.2 kb). From the GWAS genotypes, 17 markers were analyzed. Only the 3′ UTR marker rs1042778 overlapped with the markers genotyped in-house. There was >98% concordance of the 1,980 samples genotyped for rs1042778 between the in-house TaqMan ™ and GWAS genotyping platforms.
Definition of linkage disequilibrium (LD) blocks Haploview (version 3.2) was used to assess HWE and to define Shown are the number of individuals with complete scores on the indicated instrument and the number of overlapping individuals with scores on multiple instruments. The ADI-R scores are provided as total and stratified by the number of verbal and non-verbal individuals LD blocks (Fig. 1). A single individual with ASD was randomly selected from each family using a random number generator to build trios suitable for Haploview analysis (Campbell et al. 2006;Campbell et al. 2008).
Analysis of association with ASD diagnosis Family-based single marker and haplotype association analyses were performed using the family-based association test (FBAT; (Horvath et al. 2001)) and haplotype-based association test (HBAT; (Horvath et al. 2004); FBAT version 1.7.2). We analyzed only empirically determined genotypes. All FBAT and HBAT analyses were performed using the empirical variance ("-e" option) because linkage has been reported for the chromosomal region containing these genes and because the empirical variance provides a more conservative estimate of association. FBAT analysis was performed with both the additive model and the recessive model because significant association has been reported with each model (Wu et al. 2005;Yrigollen et al. 2008). HBAT was performed only on defined LD blocks, only with the additive model, and with minimum haplotype frequency set to 0.01.
Analysis of association with narrow autism diagnosis On February 23, 2010, the AGRE pedigree file was down-loaded from the AGRE website (www.agre.org). We identified 975 individuals from 618 independent families with narrow autism, defined by clinical diagnosis confirmed by both ADI-R and ADOS classification. FBAT and HBAT analyses were repeated, again using the empirical variance ("-e") option, on this smaller sample of individuals with narrowly defined autism.
Concordance and correlation among phenotype scores On February 23, 2010, six files were downloaded from the AGRE website: the phenotype scores on the SRS, the ADI-R, and each of the four ADOS modules. The total number of individuals with phenotype scores for each instrument and the number of individuals with scores on multiple instruments are present in Table 2. The phenotype scores downloaded from the AGRE website were converted directly to phenotype (.phe) files used for FBAT software analysis. Scores reported here are: (a) quantitative summation scores from individual items on the ADI-R and ADOS, (b) binary cutoff scores from the ADI-R and ADOS, (c) factor scores from a previously published principal components analysis (PCA) of the AGRE sample (Frazier et al. 2008), and (d) T scores derived from total score and subscales of the SRS. For individuals with SRS scores from  For each marker, the location on chromosome 3 and the distance from the transcription start (TrxSt) site of OXTR are indicated. Also listed is whether the marker was genotyped by the Levitt lab in-house using Taqman assays or by the Broad Institute using the Affymetrix 5.0 GWAS platform both parent and teacher scales, the SRS T scores were averaged from all available informants (Constantino et al. 2009). We performed PCA of the AGRE ADI-R scores (Campbell et al. 2010) and obtained factor structures that were indistinguishable from those previously reported (Frazier et al. 2008); therefore, we also report association analysis of two-factor scores on the ADI-R. Frazier et al. (2008) report a two-factor solution with high loadings of seven variables on the first factor (SOC1T_CS, SOC2T_CS, SOC3T_CS, SOC4T_CS, COM1T_CS, COM4T_CS, and COM2VTCS) and three variables on the second factor (COM3VTCS, BEH1T_CS, and BEH2T_CS). Four distinct modules of the ADOS are administered, depending upon age and verbal abilities of the subject, and it is not appropriate to collapse quantitative phenotype scores across the four ADOS modules. Further, there was not sufficient power in this sample to compare across ADOS modules as the number of subjects in each group ranged from 101 to 540 (mean± standard deviation=347±192; data not shown). Therefore, we did not analyze quantitative traits on each of the ADOS modules. However, each ADOS module contains a module-specific algorithm for determining whether an individual meets criteria for autism or ASD and reports a binary cutoff (yes or no) for meeting the criteria for diagnosis. The ADI-R and ADOS provide categorical variables that define thresholds for ASD diagnosis. The concordance among these binary categorical variables was determined by calculating the number of genotyped individuals with matching categorical scores compared to the total number of individuals with genotypes. The ADI-R and SRS provide continuous total scores for each individual that approximate normal distributions. Correlations among the continuous variables were calculated using StatView software.
Analysis of association with phenotype scores FBAT analysis of the phenotype scores was performed using the empirical variance ("-e" option) with no offset. Phenotype analysis was restricted to three markers, and all analyses performed are reported in the "Results" section. FBAT analysis was used for ADI-R total scores, Frazier factor combined scores, and SRS total scores, all of which approach normal distributions in this sample. FBAT was also used to analyze cutoff scores on the ADI-R and ADOS, all of which are binary (yes/no). We did not attempt analysis of individual questions on the three phenotypic instruments as the scores on individual items often fail to be normally distributed. The reported two-sided P values represent only those with positive association of the indicated allele; the indicated allele is positively associated with ASD diagnosis or phenotype.
Corrections for multiple comparisons We report in the text uncorrected P values. These results should therefore be interpreted with great caution. None of the reported associations would survive Bonferroni correction for the 25 SNPs analyzed in this study. We note, however, that this is an attempt to find associations consistent with previously reported associations of common OXTR markers with ASD risk.

Results
Linkage disequilibrium (LD) structure Analysis of LD structure indicated five LD blocks in the 39 kb that includes the OXTR gene and 10 kb of flanking sequence (Fig. 1). One LD block contains considerable overlap with the neighboring CAV3 gene, which maps 3.5 kb downstream of OXTR in a tail-to-tail arrangement. Interpretation of markers with LD overlap between the OXTR and CAV3 genes must include the possibility that genetic variants in CAV3 may contribute to ASD risk.
Association of OXTR markers with ASD diagnosis In the entire data set of 25 markers genotyped in 1,238 families including 2,333 individuals with ASD, two markers had significant association with ASD diagnosis (Fig. 2). First, the 3′ UTR marker rs1042778 G allele was associated with ASD risk using the FBAT recessive model (P=0.037). Association of the rs1042778 G allele is consistent with two previous reports of positive association (Lerer et al. 2008;Jacob et al. 2007) and further implicates a putative functional variant near the 3′ UTR of OXTR in ASD risk. Second, the rs7632287 G allele was associated with ASD using both the recessive (P=0.016) and additive (P=0.016) models, consistent with the recent report of Tansey et al. (2010). The rs7632287 marker is intergenic between OXTR and the neighboring CAV3 gene and shares LD with markers in both genes. Therefore, association of rs7632287 suggests the presence of a functional variant that affects OXTR and/or CAV3. HBAT did not reveal significant global association of any of the five LD blocks.
Association of OXTR markers with narrow autism diagnosis Family-based analysis of the 975 AGRE individuals with narrow autism (clinical diagnosis confirmed by both ADI-R and ADOS) revealed association of three markers (Fig. 3). Both variants that were associated with ASD also were associated with narrow autism diagnosis. The 3′ UTR marker rs1042778 G allele was associated with narrow autism using the recessive model (P=0.041) and the intergenic rs7632287 G allele was associated with narrow autism using both the additive model (P=0.031) and the recessive model (P=0.004). In addition, the OXTR intron 3 marker rs2268493 T allele was associated with narrow autism diagnosis using the recessive model (P=0.043), consistent with a previous positive association (Yrigollen et al. 2008). HBAT analysis did not reveal significant global association with risk for narrowly defined autism.    Fig. 1. The OXTR intron 3 marker rs2268493 (marker 12) T allele was associated with narrow autism diagnosis using the recessive model (P=0.043). The 3′ UTR marker rs1042778 G allele (marker 20) was associated with narrow autism using the recessive model (P=0.041). The intergenic marker rs7632287 (marker 23) G allele was associated with narrow autism diagnosis using both the recessive model (P=0.004) and the additive model (P=0.031) Analysis of phenotype data The ADI-R, ADOS, and SRS provide two distinct types of variables. First, the ADI-R and ADOS provide categorical variables that define binary (yes/no) thresholds for diagnosis or a phenotype domain cutoff. In our sample, the concordance (rate of agreement on an individual eclipsing the threshold) among the 15 categorical variables was high, ranging from 0.732 to 1.000 (Table 4). High concordance rates among these categorical variables should be expected for a sample selected for an ASD diagnosis. Second, the ADI-R and SRS provide continuous total scores that give a more complete description of the phenotype domains. In our sample, the correlation among these continuous variables ranges from 0.012 to 0.983, with some correlations being negative ( Table 5). The broad range of correlations among the continuous variables reflects the phenotypic heterogeneity of ASD.
Analysis of phenotype data The number of individuals with available phenotype data with each of the three instruments for measurement of ASD traits is detailed in Table 2. Phenotype analysis was restricted to the three OXTR markers that showed significant association with ASD and narrow autism risk: intron 3 marker rs2268493, 3′ UTR marker rs1042778, and intergenic marker rs7632287. Table 6 describes association of these OXTR markers with categorical variables measured on the ADI-R and ADOS. Table 7 describes association of the three OXTR markers with continuous variables on the ADI-R and SRS. Association of a genetic marker with two highly correlated phenotype measures may reflect an inability to dissect the phenotypes, rather than a contribution of the genetic variant to each phenotype independently.
Association of OXTR markers with categorical cutoff scores on the ADI-R and ADOS The intron 3 marker rs2268493 T allele was associated with 14 of the 15 categorical cutoff scores on the ADI-R and ADOS using the FBAT recessive model (Table 6). Using the FBAT additive model, the intron 3 marker rs2268493 T allele was associated only with the non-verbal communication cutoff (Table 6). The 3′ UTR marker rs1042778 G allele was not associated with any cutoff on the ADI-R (Table 6), but was associated with autism diagnosis (P=0.013), ASD diagnosis (P=0.043), autism social cutoff (P=0.036), and the autism communication plus social cutoff (P=0.021) on the ADOS using the FBAT recessive model. The intergenic marker rs7632287 G allele was associated with 14 of 15 categorical cutoff scores on the ADI-R and ADOS (Table 6). In contrast to the intron 3 marker rs2268493 T allele, the association of the intergenic marker rs7632287 G allele was independent of FBAT model used ( Table 6). The only categorical variable with which the rs7632287 G allele was not significantly associated was the non-verbal communication cutoff (Table 6), the phenotype measure with the least power in our sample ( Table 2).

Association of OXTR markers with continuous scores on the ADI-R and SRS
The intron 3 marker rs2268493 T allele was associated with the social total score (P=0.011), the non-verbal communication total score (P= 0.026), the development total score (P=0.022), and Frazier Factor 1 (P=0.011) using the FBAT recessive model (Table 7). The 3′ UTR marker rs1042778 G allele was not associated with any continuous variable score (Table 7). The rs7632287 G allele was associated with 12 of 14 continuous phenotype scores ( Table 7). The rs7632287 G allele was associated with the social total (P=0.022), the verbal communication total (P=0.010), the behavior total (P=0.002), the development total (P=0.029), Frazier Factor 1 (P=0.018), and Frazier Factor 2 (P=0.001) on the ADI-R using the FBAT recessive model (Table 7). Consistent with the categorical cutoff scores for the ADI-R, the rs7632287 G allele was not associated with the non-verbal communication total or the communication total scores on the ADI-R (Table 7) but was associated with the other six variables independent of the FBAT model used (Table 7). The rs7632287 G allele was associated with the SRS total score (P=0.029) and each of the SRS sub-scale scores using the FBAT recessive model (Table 7).

Discussion
The present study represents the largest family-based analysis to date of association of common OXTR genetic variants with ASD risk. As in previous reports of genetic association for OXTR, we found evidence of association in two distinct regions of the gene: intron 3 and the 3′ UTR. Associations of markers in intron 3 implicate the OXTR gene specifically in ASD susceptibility. The genetic associations in the 3′ UTR region, however, could influence either OXTR or CAV3 because the LD block overlaps regions of both genes. Our results are consistent with (Sullivan 2007) previous reports of association of the intron 3 marker rs2268493, the 3′ UTR marker rs1042778, and the intergenic marker rs7632287. These data also indicate, for the first time, association of OXTR genetic polymorphisms with social aspects of ASD. None of the reported nominally significant associations would survive appropriate correction for multiple comparisons, so these data should be interpreted cautiously. However, the additional suggestive evidence in support of previous genetic association data from smaller, independent family cohorts, as well as the biological plausibility, suggest that polymorphisms of OXTR may contribute to ASD risk in subsets of families that may exhibit other unique features  (Sullivan 2007;Lucht et al. 2009). The additional genetic evidence highlights two specific regions, intron 3 and the 3′ UTR, of OXTR that warrant further research.
Three recent genome-wide association studies have each identified genome-wide significant association signals at single loci on chromosomes 5p14.1 ), 5p15 (Weiss et al. 2009), and20p13 (Anney et al. 2010).  Similarly, although genome-wide copy number variation analyses have not highlighted OXTR ), a recent report described deletion of the region including OXTR, CAV3, and three neighboring genes in a single proband in one of 119 families (Gregory et al. 2009). The same report also described altered patterns of methylation in the OXTR promoter in individuals with ASD and a decreased expression of OXTR in postmortem brains of individuals with ASD (Gregory et al. 2009). Therefore, there may be multiple modes of disrupting OXTR that result in decreased expression of the oxytocin receptor and an increased risk for ASD. A bioinformatics analysis indicated that the G allele of rs7632287 may alter transcription factor binding. Consite (http://asp.ii.uib.no:8090/cgi-bin/CONSITE/consite/) predicts that COUP-TF is the only transcription factor that will bind the A allele but that the G allele will also be bound by the N-MYC, ARNT, and USF transcription factors. Experimental evidence will be necessary to determine if any of these transcription factors alter expression. A meta-analysis of OXTR association in all published data will provide a definitive answer to the association of common variants in OXTR with ASD. Our data set is publicly available at www.agre.org.
Our data suggest that the rs7632287 variant contributes to multiple domains of ASD. It should be expected that a genetic marker associated with one phenotype score will also be associated with a second highly correlated phenotype score. However, the rs7632287 G allele was associated with both Frazier Factor 1 and Frazier Factor 2 scores, even though these two phenotypic variables are only 0.070 correlated. Similarly, the rs7632287 G allele was associated with both ADI-R development total score and the SRS total score, despite the two variables being only 0.199 correlated. These data suggest that, rather than association due to highly correlated scores, the rs7632287 G allele is associated with each phenotypic domain. We were somewhat surprised by the data indicating that all domains of autism were correlated with genotype, rather than the expected statistical enrichment for OXTR relationship with specific social endophenotypes. One possible explanation for this broad association may lie in the robust developmental expression of OXTR in many more brain areas than in the adult (Shapiro and Insel 1989;Snijdewint et al. 1989;Wang et al. 1997). Thus, there could be an early role for OXT in developing neural circuitry that underlies behaviors beyond those of social interactions. As has been reported for many genes, early perturbation can lead to later onset brain dysfunctions that would not be reflected in adult expression patterns (Thompson and Levitt 2010). This also suggests that the broader functions influenced through developmental mechanisms would not be captured by intranasal OXT in healthy adults, in which functions associated with social domains show improvement. Relevant to the present study, the more widespread developmental expression of OXTR may contribute to broad domains of social phenotypes, underlying the broad association of OXTR genotypes with behavioral traits in ASD. Our data also suggest that clinical trials with OXT in young children should examine functional domains beyond those related to social behavior.