Exome sequencing of 85 Williams–Beuren syndrome cases rules out coding variation as a major contributor to remaining variance in social behavior

Abstract Background Large, multigenic deletions at chromosome 7q11.23 result in a highly penetrant constellation of physical and behavioral symptoms known as Williams–Beuren syndrome (WS). Of particular interest is the unusual social‐cognitive profile evidenced by deficits in social cognition and communication reminiscent of autism spectrum disorders (ASD) that are juxtaposed with normal or even relatively enhanced social motivation. Interestingly, duplications in the same region also result in ASD‐like phenotypes as well as social phobias. Thus, the region clearly regulates human social motivation and behavior, yet the relevant gene(s) have not been definitively identified. Method Here, we deeply phenotyped 85 individuals with WS and used exome sequencing to analyze common and rare variation for association with the remaining variance in social behavior as assessed by the Social Responsiveness Scale. Results We replicated the previously reported unusual juxtaposition of behavioral symptoms in this new patient collection, but we did not find any new alleles of large effect in the targeted analysis of the remaining copy of genes in the Williams syndrome critical region. However, we report on two nominally significant SNPs in two genes that have been implicated in the cognitive and social phenotypes of Williams syndrome, BAZ1B and GTF2IRD1. Secondary discovery driven explorations focusing on known ASD genes and an exome wide scan do not highlight any variants of a large effect. Conclusions Whole exome sequencing of 85 individuals with WS did not support the hypothesis that there are variants of large effect within the remaining Williams syndrome critical region that contribute to the social phenotype. This deeply phenotyped and genotyped patient cohort with a defined mutation provides the opportunity for similar analyses focusing on noncoding variation and/or other phenotypic domains.

Understanding contributions to social phenotypes in particular for WS may define genes that regulate human social behavior, providing insight not only into WS, but also in other disorders as well as possible modifiers of social behavior in the general population. Deleting one copy of the genes in the WSCR produces the personality profile observed in WS, which consists of prosocial behaviors such as gregariousness, empathy, retained expressive language skills, and low levels of social anxiety, in spite of high anxiety in other domains (Doyle, Bellugi, Korenberg, & Graham, 2004;Gosch & Pankau, 1997;Järvinen et al., Dup7 (Klein-Tasman & Mervis, 2018;Richards, Jones, Groves, Moss, & Oliver, 2015). The similarities and differences in the social communication domains of WS and ASD have been described, and suggest that while both disorders show deficits in social communication, the WS group was not as impaired as the ASD group (Klein-Tasman et al., 2007. Unlike ASD, there is no sex bias in the frequency of WS and severity of social and cognitive phenotypes are similar across both sexes (Brawn & Porter, 2014;Dykens, 2003).
As in many diseases of haploinsufficiency, within WS there remains considerable variability in expressivity of the phenotypes, despite the very homogeneous genetic cause. It is thought that both genetic background and the environment introduce variation in the expression of a phenotype. The fact that individuals with WS are hemizygous for 26-28 genes has led to the assertion that variation in the remaining allele could contribute to the severity of symptoms in WS (Delio et al., 2013;Merla, Brunetti-Pierri, Micale, & Fusco, 2010). The presence of only one copy of genes in the WSCR could unmask the effects of recessive alleles in the region that are more difficult to detect in a diploid setting. Indeed, this logic has been applied to investigate the variability in the cardiovascular phenotype. Delio et al., 2013 sequenced the exons that make up the ELN gene in a sample of 55 individuals with WS, but found no clear link between severity of phenotype and remaining genetic variation. However, no similar studies have investigated the social profile of WS, in spite of the fact that there is some evidence that common variation in the region can influence social behavior in the general population. For example, variation in the GTF2I gene has been associated with the WS cognitive profile, autism, oxytocin reactivity, amygdala activity, and social anxiety (Crespi & Hurd, 2014;Jabbi et al., 2015;Procyshyn, Spence, Read, Watson, & Crespi, 2017). Furthermore, genes outside of the WSCR are also likely to affect aspects of social behavior. In particular genes that are associated with ASD have a profound effect on social interaction and could harbor variants that modify the phenotype of individuals with WS.
Here, we employ whole exome sequencing to understand how genetic variation within the WSCR, and other protein coding genes, impacts the severity of the WS social phenotype. We generate a rich catalog of genetic variants identified from 85 individuals with the typical WS deletions; each individual has also been assessed with the Social Responsiveness Scale-2 (SRS) questionnaire, a quantitative measure of reciprocal social behavior. The SRS was first developed to quantify autistic traits in both the general and clinical populations Moreno-De-Luca et al., 2015). SRS scores have also been used to describe different aspects of the social phenotype in WS (Klein-Tasman et al., 2010). We then employ a three-tiered approach to screen for the existence of alleles that contribute to SRS scores in the context of a potentially sensitizing WSCR deletion, ordering the analyses to conserve statistical power. First, we describe the genetic variants observed in the remaining WSCR and test if they can explain the variance in the SRS scores. We find little evidence that these common or rare variants in the region are associated with SRS scores. Next, we go beyond the WSCR and test variants in 71 genes known to be associated with ASD (Sanders et al., 2015), reasoning variation that contributes to autistic features in non-WS children may modify autistic features in the WS cohort as well. Finally, we test variants throughout the whole exome. We find no genetic variants of sufficient effect size to support the hypothesis that they contribute to the social phenotype in this sample of individuals with WS. However, we have more thoroughly described the variation in the WSCR region as it relates to social behavior and provide the largest genetic dataset to date of individuals with typical WS deletions for future analyses of other phenotypic domains.

| Ethical compliance and samples
This study was conducted with approval of the IRBs at Washington University School of Medicine and the National Institutes of Health. Consent was obtained prior to inclusion in the study. Once enrolled, participants provided a DNA sample by blood or saliva and their care-givers filled out health related questionnaires. The 85 individuals that make up our sample have ages that range from 2.5 to 65.5 years with a mean of 16.1 years. Caregivers provided a self-reported ethnicity. The majority of the sample was reported as white (77 individuals). There are two individuals that are African American, three Chinese, and three others.

| Confirmation of diagnosis
WS diagnosis and typical deletion size was confirmed using either chromosomal microarray or quantitative PCR. In some cases, clinical microarray results were derived from the medical record. Array type varied by individual. For the remaining individuals, some received a research array (Cytoscan HD, Applied Biosystems) with analysis using the accompanying ChAS software. Others underwent deletion size assessment using quantitative PCR for genes within and outside of the Williams region using Taqman copy number probes (Thermo-Fisher, AUTS2: Hs04984177_cn, CALN1: Hs04946916_cn, FZD9: Hs03649975_cn, CLIP2: Hs00899 301_cn, HIP1: Hs00052426_cn, POM121C: Hs075298 20_cn). Copy number analysis was done according to the manufacturer's instructions and output data analyzed using their Copy Caller software. All individuals were confirmed to have deletions that included the WSCR genes ELN, FZD9 and CLIP2, but did not include genes external to the typical deletion such as CALN, AUTS2, POM121C or HIP1 (data not shown).

| Social responsiveness scale
The social responsiveness scale-2 (SRS) is a 65-item questionnaire that measures aspects of social interaction that make up the core symptoms of autism spectrum disorders. The output is a total raw score as well as a T-score that is adjusted for sex, age, and the relationship of the reporter to the proband. The total score is made up of the scores of five subcategories that are impaired in ASDs: social awareness (AWR), social cognition (COG), social motivation (MOT), social communication (COM), and behaviors typical of autism such as restricted interests and repetitive behaviors (RRB). The response to each question ranges from 1 (not true) to 4 (almost always true). The T-scores are binned into four groups: normal <59, mild between 60 and 65, moderate between 66 and 75, and severe >76. For this study, the agespecific (preschool, school age, or adult) SRS-2 was completed by the participant's caregiver and analyzed as a Tscore that is adjusted for sex, age, and the relationship of the reporter. We provide values from the general population that have been previously reported for comparison .

| Sequencing and variant calling
Whole exome sequencing and alignment was performed at Washington University in St. Louis by the McDonnell Genome Institute on 85 DNA samples from individuals with WS. Exomes were captured using Nimblegen SeqCap EZ Choice HGSC Library version 2.1, which targets 45.1 Mbp covering 23,585 genes and 189,028 nonoverlapping exons. Exomes were aligned to the GRCh37-lite genome using bwa -mem v0.7.10 ) default settings, samtools v0.1.19 ) was used to assign mate pairings, sort, and index the bam files. Duplicates were marked using Picard MarkDuplicates v1.113.
Variant calling was done following GATK best practices on the aligned exomes . Briefly, using GATK v3.6.0 indels were realigned and the base quality scores recalibrated. Variants were initially called per sample using the haplotype caller tool, followed by jointly calling variants. To improve variant calls, we recalibrated variants and used a truth sensitivity tranche of 97 for SNPs, and a truth sensitivity tranche of 94 for indels. These thresholds were chosen to maximize the number of known and novel variants while still being stringent enough to limit the number of false positive variant calls. To further filter the variants we used the VariantFiltration tool to filter variant sites that had lower than a 10× average coverage or an inbreeding coefficient < −0.20 to remove sites with excess heterozygosity. Genotype calls were filtered and considered to be missing if they had a genotype quality score of <20, which refers to a 99% probability that the call is correct. Finally, using vcftools v0.1.14 (Danecek et al., 2011), we removed sites that had a genotype missing rate of >10%, as well as sites that no longer showed any variation. This produced a call set of 202,820 variant sites. The final call set has a Ti/Tv ratio of 2.76 and a dbSNP rate of 88.5%. These metrics are consistent with quality variant calls and a low false positive rate.

| Variant annotation
The variant call set was split into three groups using vcftools: (1) variants in the Williams syndrome critical region (WSCR) defined by hg19 coordinates chr7:72,395,660-74,267,841 (2) variants located in 71 genes associated with ASD (Sanders et al., 2015), and (3) the remaining nonoverlapping variants. All sets include exonic variants as well as variants located in introns that are pulled down by the capture reagents. Bcftools v1.2 ) was used to split multiallelic sites into separate lines for each allele and left normalized so positions would be compatible with ANNOVAR annotation files version 2016-02-01 (Wang, Li, & Hakonarson, 2010). The ANNO-VAR table_annovar.pl function was used to annotate all three variant call sets with the RefSeq gene annotation, variant consequence, ExAC allele frequency (Lek et al., 2016), sample specific allele frequency, dbsnp147 name, clinical significance assessed by ClinVar (Landrum et al., 2016). Missense variants were also annotated with measures of deleteriousness compiled in dbNSFPv3.3a (Liu, Wu, Li, & Boerwinkle, 2016). We highlight the CADD PHRED score and MetaLR as two measures of variant deleteriousness. CADD scores are defined at each base in the genome and for every possible single-nucleotide change (Kircher et al., 2014). CADD scores compare 65 annotations, including functional data as well as conservation scores, between fixed human derived alleles and simulated variants. Deleterious variants should be depleted in the observed fixed alleles and not in the simulated variants. CADD PHRED scores represent the relative rank of a CADD score compared to all other possible allele CADD scores; a CADD score of 10 means this allele is ranked as the top 10% of all possible CADD scores. Larger CADD PHRED score indicates an increased predication of deleteriousness. MetaLR uses logistic regression to incorporate information from nine other variant annotations that consider function as well as conservation (Dong et al., 2015). The model was trained on true deleterious variants and true neutral variants described in the Uniprot database. The composite MetaLR score was found to have greater predictive ability than any of the single scores that make up MetaLR.

| Power analysis
We performed a power analysis to provide the limits of genetic effects that we would be able to detect given our cohort size. For future studies we also calculate the sample sizes that would be needed to detect different magnitudes of genetic effects. We used the Genetic Power Calculator (Purcell, Cherny, & Sham, 2003). We calculated the predicted power of the current sample size n = 85 using a pvalue threshold corresponding to the Bonferroni corrected alpha for each set of analyses (WSCR 34 variants, alpha = 0.00147, ASD 381 variants, alpha = 0.000131, WEX 66620 variants, alpha = 7.5 × 10 −7 . Our main hypothesis is variants on the remaining WSCR allele affect the social phenotype; we wanted to calculate the sample sizes that would be required to detect different size genetic effects in the WSCR at different levels of power. We again used the alpha threshold based on the 34 common variants we identified in the exons of the WSCR and report the sample size required to achieve a specific power.

| Common variant analysis
The variant call files were converted to plink binary bed format using the GATK tool VariantToBinaryPed. We used PLINK v1.9 (Purcell et al., 2007) -linear option to conduct a quantitative trait association using the SRS T-score as the quantitative trait. Ancestry was controlled for by including the first four principle components, determined by the -pca function in PLINK, as covariates along with sex and age. We used alleles that had a minor allele frequency (MAF) of 0.05 or greater. We performed the association analyses on the three separate groups of variants described in the previous section. It should be noted that allele frequency in the Williams syndrome critical region is inflated because of the hemizygous state of the region in individuals with WS. A MAF of 0.05 in this region corresponds to an allele count of four. In all cases we report the effect size of a variant under an additive model. Though the small sample size of this study limits power, in an exploratory fashion we also performed the same quantitative trait analysis on each of the subscores of the SRS using variants in the WSCR, ASD genes, and the whole exome.

| SKAT-O
SKAT-O (Lee, Wu, & Lin, 2012)was implemented in the R v3.1.3 environment. SKAT-O fits a multiple linear regression of all SNPs located in a user provided region. The framework in SKAT-O allows for correlation between SNPs in a region, where if all SNPs are perfectly correlated this would become a burden test, but also allows SNPs in the same region to have effects in opposite directions. Significance is assessed by region rather than by SNP. We considered each gene that harbors a variant in the WSCR as a separate region for a total of 26 regions. To test for an overall effect of variants in the ASD genes we collapsed the 61 autosomal genes into one region. We used the beta function shape parameters (1,50) to put more weight on SNPs that have lower minor allele frequency, reasoning that rare causal alleles potentially have a greater effect size. We again controlled for age, sex, and the first four principal components.

| Polygenic risk score
Polygenic Risk Scores (PRS) can be used to test if there is a contribution of many loci of small effect on the phenotype of interest by summing the effects of variants that may have not reached genome-wide significance. For a discovery set, we used the publically available summary statistics from the most recent Psychiatric Genome Consortium genome-wide association study (GWAS) of autism spectrum disorder (The Autism Spectrum Disorders Working Group of the Psychiatric Genomics Consortium 2017), reasoning that genetic risk for autism would contribute to SRS scores. The best-fit PRS was determined using the high-resolution functionality in the PRSice software (Euesden, Lewis, & O' Reilly, 2015). All the variants identified throughout the exome with a MAF >0.05 and that are also present the in the discovery set were used to calculate the PRS. Sex, age, and the first four PCs were included as covariates. After clumping there were a total of 23,191 variants used to calculate the PRS. PRSice was used to calculate the significance of the PRS at the best-fit p-value threshold using 10000 permutation to determine an empirical p-value. PRS for each of the samples was calculated for the total SRS T-score as well as the subscores.

| Other statistical analyses
All remaining statistical tests were done in the R v3.1.3 environment. Two sample t-tests were used to compare the means of two groups. ANOVA was used to test differences in mean of subscales of SRS. TukeyHSD post hoc comparison was performed using the multcomp package. The qqman (Turner, 2014) package was used to generate manhattan and qq plots.
We examined the SRS and its subscores in depth. In our sample, the SRS T-scores are continuously distributed in the WS population with a male mean T-score ± SD of 64.58 ± 12.28 (mean male raw score ± SD 74.53 ± 32.03) and female mean T-score ± SD of 62.94 ± 11.04 (mean female raw score ± SD 67.08 ± 26.04) (Figure 1). There is no significant difference in SRS T-scores (t 70.76 = 0.6365, p = 0.52) or raw scores (t 65.907 = 1.1445, p = 0.257) between sexes. To benchmark the WS values, Constantino & Todd, 2003 measured raw SRS scores in 788 twin pairs from the general population ranging in ages between 7 and 15 and estimated the mean male raw score ± SD as 35.3 ± 22.0 and the female mean raw score ± SD as 27.5 ± 18.4; males and females were significantly different. In our analysis, we show that individuals with WS have SRS scores that are shifted toward the more impaired end of the spectrum, and we do not detect any significant sex differences in WS, which has been observed in the general population.
Our results largely replicate the results seen in Klein-Tasman et al., 2010;. The overall T-score distribution reveals that 40% of our samples fall into the no clinically significant impairment range, followed by 41.1% with mild to moderate deficits, and 18.9% with severe deficits. The number of individuals showing no clinical signs in our sample is higher than the 13.4% observed when the parents completed the SRS in Klein-Tasman et al., 2010; but more similar to the teacher reported results of 38.8%. The subscores also follow a similar pattern to what has been reported previously (Klein-Tasman et al., 2010). There is a significant effect of subscale on the T-scores (F 4,420 = 24.759, p < 0.001) (Figure 1b). Post hoc Tukey all-pairwise comparisons show that social motivation has significantly better T-scores than all other subscales, consistent with Klein-Tasman et al., 2010. The social awareness and communication scales are not different from each other, but both show less impairment than social cognition and restricted and repetitive behaviors. Social cognition and restricted and repetitive behaviors were significantly more impaired than all other subscales, but not each other.
The distribution of SRS scores in WS points to the possibility of additional genetic variants that modify the social phenotype. First, we see a larger standard deviation in the SRS data in our sample compared to that of the norming population from . The extra variance suggests individuals with WS are more sensitive to genetic or environmental factors that modify social behavior. Second, in our sample there are only two individuals that show severe social motivation deficits, and these individuals also show severe deficits in the total SRS T-score as well as all other subscales. These outliers also suggest some individuals may harbor additional rare variants of large effect size resulting in a phenotype that is more frankly autistic. To test these two hypotheses, we generated and analyzed exome sequences from this cohort of WS patients.

Williams syndrome critical region
Williams syndrome individuals are hemizygous for 1.5-1.8 Mbp on chromosome 7q11.23. Since they only have one remaining allele, our primary hypothesis was that second hits in genes believed to impact social phenotypes within the WSCR would produce more extreme social phenotypes. We performed whole exome sequencing on 85 individuals, all whom have an SRS score. We called 120 variants in the remaining WSCR and annotated them with the allele frequency in our sample, ExAC allele frequency, mutation consequence, clinical significance as assessed by ClinVar, and scores that assess deleteriousness of missense variants cataloged in dbNSFP. (Supporting Information  Table S1). Table 1 shows the 55 exonic variants discovered in the region. For display purposes we have only included the CADD PHRED score and the MetaLR score, which is a composite score that incorporates information from nine other measures of deleteriousness and has been shown to have more predictive power than the individual component scores (Dong et al., 2015).
We first examined this set of variants to determine if any loss-of-function variants might be present in individuals with particularly severe SRS scores in our sample. Upon inspection of the exonic variants, we notice no severe likely protein-truncating variants. As homozygous nulls for at least two genes in this region (ELN and GTF2I) are expected to be lethal (Li et al., 1998;Sakurai et al., 2011), we also assessed missense mutations in these genes that might alter function. Based upon the predictions of MetaLR all the missense mutations called are expected to be tolerated. None of the variants were reported as pathogenic in ClinVar. The highest CADD scores observed are a Severity bins of SRS and subcategory scores novel variant and SNP rs35607697, both located in the TBL2 gene. Another novel variant was identified as a synonymous change in the BAZ1B gene. Similar results are found for nonexonic variants in the region (Supporting  Information Table S1). This suggests that beyond the reduced copy number of the entire WSCR, neither a second rare deleterious coding variant nor any common missense mutations in the WSCR explain individuals with outlier SRS scores. It should be noted that we did not identify any variants in GTF2I, one of the primary candidates for mediating the social-cognitive profile.

| Association analyses
To test the hypothesis that individual variants in the WSCR can explain the variance in the SRS scores in our sample, we perform classic quantitative trait loci associations. Rare disease populations by definition will have small sample sizes such as in this study. We calculated the power of our study to be able to detect variants with different effect sizes and also calculated the number of samples that would be needed to reach a certain power given an effect size (Figure 2). We calculated the power for analyzing variants in the WSCR, variants in 71 ASD genes, and the remaining variants identified throughout the exome. Since we are conducting fewer tests in the WSCR, we have the most power in this analysis, however we are still only powered to detect very large effect sizes that might be unmasked by the hemizygosity of the region, such variants would need to explain more than 10% of the heritability of the trait to achieve 80% power. Most effect sizes for common variants in diploid regions of the genome typically assessed by GWAS for complex traits explain around 1% of the heritability of the trait (Manolio et al., 2009). In order to be able to detect variants that explain 5% of the variance of the trait with 80% power using only variants in the WSCR would require 312 individuals (Figure 2b). We then performed a quantitative trait association analysis of common variants in the WSCR on the SRS T-scores from the whole cohort. We used PLINK to test for an association on each of the 34 common variants in the WSCR, defined as MAF > 0.05, which corresponds to an allele count of at least four in the WSCR due to the hemizygosity of the region. We adjusted for age, sex, and ancestry. We found no association between any SNP and SRS that survived multiple comparison corrections (Figure 3a). The top five SNPs are displayed in Table 2. Interestingly, the most significant SNP, rs2074754, is located in the BAZ1B gene, which has been previously implicated in contributing to the cognitive phenotypes in WS (Lalli et al., 2016). Furthermore, the next most nominally significant SNP is rs61438591, an intronic variant in the GTIF2RD1 gene, another gene highly implicated in the cognitive and social phenotypes seen in WS (van Hagen et al., 2007;Howard et al., 2012;Schneider et al., 2012;Young et al., 2008).
As the common variants in WSCR showed no association, we wanted to test for the possibility that rare variants could contribute to the variability in SRS T-scores. To test this, we used SKAT-O, which tests all variants in the region at once and weights each variant by its minor allele frequency. Similarly, we included age, sex, and ancestry as covariates. We tested each gene in the WSCR independently, because we hypothesized only certain genes in the region, such as STX1A, LIMK1, CYLN2, BAZ1B, GTF2IRD1 (Dai et al., 2009;Fujiwara et al., 2016;Gao et al., 2010;van Hagen et al., 2007;Hoogenraad et al., 2002;Lalli et al., 2016;Meng et al., 2002;Morris et al., 2003;Sakurai et al., 2011) that have been implicated in the cognitive phenotypes would contribute to the social phenotype rather than the entire region. While no gene p-value survives multiple testing corrections, the ELN gene has the most nominally significant p-value of 0.013.
The results of our analysis of variation in the WSCR suggest that common and rare variants in the remaining allele do not strongly influence social behavior in WS. This does not exclude the possibility that a second deleterious hit or common variation in other genes outside the region contributes to the variation in the SRS T-scores. To test this, we next examined variation in 71 genes known to be associated with autism spectrum disorders (Sanders et al., 2015). These genes should be enriched for loci that affect social behavior and genetic variation in these genes could contribute to variability seen in WS. We called 1,367 variants in the 71 genes (Supporting Information Table S2). We annotated the variants as above, with clinical significance and measures of deleteriousness compiled in dbNSFP. There are 313 (22.9%) variants that had at least one submission to ClinVar. None of these variants had previous evidence to support pathogenicity. There are 33 missense variants predicted to be deleterious by MetaLR that are seen in 36 individuals in our sample. Despite having a putatively deleterious variant the distribution of SRS Tscores is similar between individuals either carrying or lacking deleterious variants in these genes (t 82.999 = 0.6878, p-value = 0.4935). There are seven variants that should result in a truncated protein, one stop gain in the USP45 gene and six frameshift mutations. Only one sample harboring one of these mutations has a severe SRS T-score of 77. All these protein-truncating mutations are also observed in the ExAC cohort.
We next tested for associations of each of the 381 common variants (MAF >0.05) in these genes. No SNP was significant after multiple testing corrections (Figure 3b). The top five SNPs are located in Table 2. Since each of these genes has been associated with ASD, we hypothesized that rare and common variants in each of the genes could contribute to SRS. We performed SKAT-O on the variants located in the autosomal ASD genes altogether, which also showed that there is little evidence to support variants in these 68 ASD genes have a strong effect on SRS T-scores, p = 0.431. While it would be underpowered for any but the largest effect sizes (Figure 2a), for thoroughness we did an unbiased scan of the whole exome. We also examined the polygenic contribution of common variants to the SRS. The common variant analysis was performed on 66,620 variants ( Figure 3c). The most nominally significant single SNP is rs527221 located in the DMPK gene, which is responsible for causing type 1 myotonic dystrophy (Brook et al., 1992) ( Table 2). While there is suggestive evidence for single variants such as rs527221, we calculated the polygenic risk scores (PRS) for each of the individuals in our sample to test if exome wide there are many SNPs of small effect that contribute to the social phenotype in WS. We used the summary statistics from the most recent PGC GWAS on autism spectrum disorders to calculate the PRS for our sample (The Autism Spectrum Disorders Working Group of the Psychiatric Genomics Consortium 2017). We reasoned the polygenic risk of autism would be correlated with the SRS because this is a questionnaire used to assess Notes. a "." Refers to information that is not applicable; b "T" the missense mutation is predicted to be Tolerated. behaviors that are affected by autism. Variants from the PGC GWAS were included if the p-value for the variant was under the threshold determined by the high-resolution screen in the PRSice software (Euesden et al., 2015). Interestingly, only the PRS for the motivation subscore was nominally significant (p = 0.033), but after permutation to determine an empirical p-value it was not significant (p = 0.308). The correlations of the PRS for each of the samples and the subscore as well as total SRS are shown in Supporting Information Fig. S1. Counterintuitively, there is a negative correlation between the PRS and motivation subscore. While this is the largest correlation between the PGS and subscores it implies that more genetic risk for autism leads to a lower and less impaired social motivation Tscore. However, given the small sample size and small number of SNPs available from whole-exome sequencing compared to whole-genome genotyping we are wary of making strong conclusions from this analysis. We and others (Klein-Tasman et al., 2010) have shown that individual subscores of the SRS are affected differently by the deletion of the WSCR. Therefore, we wanted to rule out the possibility that variants are indeed affecting specific subscales of social behavior, but that testing the total SRS score is masking those effects. Thus, in an exploratory manner, we repeated the quantitative trait loci associations for each of the subscores of the SRS using the variants in the WSCR, 71 ASD genes, and the remaining whole exome variants. Since the sample size is small we conducted these associations for exploratory and hypothesis generating purposes. The top five SNPs from each association are reported in Supporting Information Tables S3-S5. For each of the analyses we see similar variants showing the highest association as were associated with the total SRS, likely due to the high correlation between the SRS and the subscores (Supporting Information Fig. S2). Thus, an analysis of the total SRS was not masking independent genetic effects on each subscale.

| DISCUSSION
Phenotypic variability has been appreciated in many of the symptom domains of WS including the cardiovascular phenotypes, the unique cognitive profile, and in social behavior (Anney et al., 2012;Joyce, Zorich, Pike, Barber, & Dennis, 1996;Porter & Coltheart, 2005). Here, we have described the variability of reciprocal social behavior in a sample of 85 individuals with the typical WS deletion using the SRS-2. Our results replicate the findings of Klein-Tasman et al., 2010, revealing that overall individuals with WS have SRS scores that are shifted to the more socially impaired end of the distribution, with most problems relating to the social cognition and restricted and repetitive behavior subscales of the SRS while social motivation is spared.
We also note that sex differences in the general population have been reported previously in the literature for SRS. These sex differences were not consistent with different genetic factors contributing to the SRS in boy and girls, but due to discrepant effects of common genetic and environmental factors on SRS, such as differences in sensitivity to environmental factors or the X-inactivation phenomenon . However, we do not see evidence of sex effects in our sample of individuals with WS. The magnitude of the difference between males and females in our sample is similar to what was reported in the general population, so our lack of a significant finding could be due to our small sample size. The standard deviation of the SRS is large in both the general population and still larger in the WS population, so it may also be that larger sample sizes are needed to overcome the considerable variance in the data. The fact that the WS population has a larger standard deviation could also suggest that individuals with the deletion are sensitized to other factors that contribute to variation in the SRS such as background genetic variation or environmental factors. We performed whole-exome sequencing on our sample of 85 individuals to test for additional genetic contributions to the variability seen in social behavior in individuals with WS. We used the identified variants to test the hypothesis that genetic variation in the remaining WSCR allele can explain some of the variability in SRS T-scores. Genes in this region have a dosage sensitive effect on social behavior evidenced from the contrasting social phenotypes of the WS deletion and the reciprocal duplication, suggesting that variants in the remaining WSCR allele that affect expression or function of the genes could further contribute to the social phenotype (Merla et al., 2010). We called 120 variants in the WSCR with 55 variants being exonic. We used evidence such as the amino acid change, clinical significance suggested by the ClinVar database, and multiple algorithms to predict the consequences of the variants. Within the WSCR we do not find any variants that cause protein truncation. None of the missense variants are predicted to be deleterious based on the MetaLR composite score. Of the nine variants that have been submitted to ClinVar, all were described as benign or likely benign. A quantitative trait association analysis using the common variants in the region resulted in no SNP that survived multiple testing corrections. The most significant SNP, rs2074754, is a synonymous SNP in the BAZ1B gene. This gene encodes for a protein product in the bromodomain protein family that modifies chromatin to affect transcription and has been implicated in the cognitive phenotypes in WS. Knocking down this gene in human derived induced pluripotent stem cells upregulates genes involved in mitosis as well as downregulating genes that are involved in the development of the nervous system (Lalli et al., 2016) The second most nominally significant SNP, rs61438591, is an intronic variant in GTF2IRD1, which encodes for a transcription factor that has been suggested to contribute to the cognitive and social behavior deficits (Dai et al., 2009;Howard et al., 2012;Schneider et al., 2012;Tassabehji et al., 2005;Young et al., 2008). If future studies with increased power replicate this association, it would suggest that noncoding variation, perhaps controlling the expression of this gene, might contribute to variation is social behavior. We also tested the association of all variants in the WSCR using SKAT-O. This test indicated no variants with sufficient effect size were detected in the WSCR.
While we have not shown evidence that variants in the remaining WSCR contribute to the social phenotype in WS, we cannot conclusively discard this hypothesis. However, our study does clearly indicate that the alleles genotyped here are either not causative or exert too small an effect size on SRS for our current power (Figure 2), but it does not rule out variants of small effect on social behavior in the region. Research on other copy number variants associated with ASDs showed that larger CNVs tended to have genes of smaller individual effect size and suggests the phenotype of the overall CNV is due to the cumulative effect of each of those genes (Sanders et al., 2015). Further we did not detect any variants in the gene GTF2I, which has been highly suspected of contributing to the social behaviors in WS (Borralleras, Sahun, Pérez-Jurado, & Campuzano, 2015;Crespi & Hurd, 2014;Dai et al., 2009;Morris et al., 2003;Sakurai et al., 2011). The lack of variant calls in our sample could be due to the fact that GTF2I is under stringent purifying selection. Indeed, looking at the ExAC data covering this gene, they show that there are fewer missense variants than expected by chance. ExAC discovered 62 synonymous and 56 missense mutations in 60,706 people (Lek et al., 2016). In our sample of 85 individuals we would expect to see variants in ExAC that have an allele frequency of >0.0059, which is an allele count of one in our sample. There are ten variants with an allele frequency >0.0059 detected in ExAC, only three of which are exonic. Thus, we would need a much larger sample size to investigate coding variants in GTF2I. The two linked variants in GTF2I that have previously been associated with oxytocin responsiveness and amygdala reactivity, rs1322743 and rs4717907, are intronic and were not covered in our sequencing Swartz et al., 2017).
We further used the genetic data to investigate the role of variation in 71 genes that have been associated with ASD. WS and ASD do show phenotypic overlap Klein-Tasman et al., 2009), and we reasoned that these genes should be enriched for functional roles in social behaviors. Likewise, the presence of outlier scores on the SRS that indicated severe impairment, suggested there could be possible second deleterious hits on top of the WS deletion in our dataset. Second hits are expected to be rare but have been observed in WS to explain a case of a child with comorbid seizures (Popp et al., 2016). Inspecting the 1,367 variants discovered in the ASD genes, 313 variants have been previously submitted to ClinVar, none of which show evidence for any pathogenicity. We observed seven protein-truncating mutations that do not associate with severe SRS T-scores. Several missense mutations were predicted to be deleterious, but there was no association between individuals that had a putative deleterious variant and a more impaired SRS score. Testing the common and rare variants in these genes showed no associations with the social phenotype. Similar results were found when we performed the association analyses on all the variants discovered in the cohort. The most significant SNP was rs527221, a nonsynonymous variant in the DMPK gene, which is responsible for causing type 1 myotonic dystrophy, severe childhood forms of which have been associated with ASD (Ekström, Hakenäs-Plate, Samuelsson, Tulinius, & Wentz, 2008). We also tested if polygenic risk for increased ASD liability is associated with the SRS T-score and subscores. This boosts our ability to detect the impact of many loci with small effects. The largest correlation was between the PRS and the social motivation subscore, although this was not significant.
WS seems to affect specific domains of social behavior as evidenced by significant differences between the subscores of the SRS. This observation led us to an exploratory examination of associations with the subscores of the SRS and test if different genetic variants contribute to each subscore. Overall using variants from the WSCR, ASD genes, or the whole exome identified the same variants as nominally significant. The SRS and the subscores are very correlated, but the social motivation in the WS sample is the least correlated with all other scores. This reflects that fact that social motivation tends to be rated within the normal range in WS, while the other scores are often higher. Interestingly, the whole exome association on the motivation T score leads to the lowest FDR values compared to the other scores, suggesting that there may be more genetic signal when using this subscale. Indeed, this decoupling of the social motivation subscale from other SRS items highlights the possibility that the social motivation subscale might provide useful clinical information going forward; individuals carrying the WSCR deletion yet not showing a spared social motivation might warrant a deeper examination for additional factors impacting their presentation.
There are several limitations to our study that should be addressed in future research. First this study genotyped and assessed only the probands and not their parents. Having genetic information from trios would allow us to distinguish between variants that are inherited or de novo, which would aid in interpretation and prioritization of variants. Furthermore, being able to compare the SRS score of the individual with WS to biparental SRS mean would let us control for effects of background genetic variation (Moreno- De-Luca et al., 2015). Second, we are limited to investigating exonic variation. While interpretation of exonic variants is more straightforward because they potentially disrupt coding sequences, and can aid in the detection of deleterious rare variants, we could be missing important regulatory information that is located in promoters or introns of genes. Third, we were not able to control for intellectual functioning of the individuals with WS. The SRS has been reported to not correlate with intellectual functioning , but Klein-Tasman et al., 2010 found significant negative correlations between intellectual functioning and the total SRS T-score when parents completed the report, but not when teachers completed the report. SRS values have been shown to be dependent on levels of expressive language, nonverbal IQ, and behavioral problems. A subset of SRS questions was selected to ameliorate these dependences (Sturm, Kuhfeld, Kasari, & McCracken, 2017). The short form of the SRS as well as other questionnaires that assess adaptive skills and social behaviors could be used in the future to provide supporting information about the social phenotype and underlying genetics in WS. Finally, while our study represents the largest single collection of WS samples reported to date, it is only powered to detect strong effects of common variants due to our small sample size. This is challenging to overcome due to the low prevalence of WS.
In conclusion, we have tested the hypothesis that variation in the remaining WSCR allele affects the social phenotype of individuals with WS, by applying whole exome sequencing to a sample of 85 individuals with typical WS deletions. We show that common and rare variants in the region do not associate with SRS T-scores in our sample. Furthermore, we show that variation outside of the region does not account for the social variability. This is not to say that genetic variation does not play a significant role in phenotypic variability in WS, but that it will require larger sample size to detect. In the future, applying whole genome sequencing to a sample of individuals with WS might elucidate the roles of genetic variation in the regulatory elements. Whole genome data could also allow for more accurate breakpoint determination. Redundant sequences in the low copy number repeat areas at either end of the WS deletion prevent accurate end point detection by CMA. This will be an interesting avenue to pursue in order to investigate how deletion size variation among individuals with typical 1.5 to 1.8 MB deletions contributes to social behavior. For example, Porter et al. showed that those with larger (1.8 Mb deletions) had decreased executive functions (Porter et al., 2012). It is also worth noting that the current genetic data set has additional clinical data available, which