Participants

The current study included participants from 23andMe (primary GWAS - SQ-R), from ALSPAC (GWAS of scores on the Social and Communication Disorders Checklist (SCDC)) and autistic individuals from the Simons Simplex Collection (SSC), the Autism Genetic Resource Exchange (AGRE), and the EU-AIMS LEAP and PARIS cohorts.

23andMe participants

Research participants in the GWAS of the SQ-R were from 23andMe and are described in detail elsewhere67,68. All participants provided informed consent and answered surveys online according to a human subjects’ research protocol, which was reviewed and approved by Ethical & Independent Review Services, an external AAHRPP-accredited private institutional review board (http://www.eandireview.com). All participants completed the online version of the SQ-R on the 23andMe participant portal. Only participants who were primarily of European ancestry (97% European Ancestry) were selected for the analysis using existing methods69. Unrelated individuals were selected using a segmental identity-by-descent algorithm70. A total of 51,564 participants completed the SQ-R (males = 26,063, and females = 25,501).

ALSPAC participants

ALSPAC is a longitudinal cohort which recruited pregnant mothers in the Avon region of the UK. The ALSPAC cohort comprises 14,541 initial pregnancies from women in Avon resulting in a total of 13,988 children who were alive at 1 year of age. Children were enrolled in additional phases, described in greater detail elsewhere71. This study received ethical approval from the ALSPAC Law-and-Ethics Committee, and the Cambridge Human Biology Research Ethics Committee. Written informed consent was obtained from parent or a responsible legal guardian for the child to participate. Assent was obtained from the child participants where possible. We conducted a GWAS of scores on the SCDC in 5,421 individuals from ALSPAC.

Other cohorts

We included data from four cohorts to conduct polygenic score and bivariate genetic correlation analyses. The SSC (n = 2221 unrelated autistic individuals) consists of simplex autistic families, and are described elsewhere72. The AGRE cohort (n = 482 unrelated autistic individuals) consists of multiplex autism families, details of which are provided elsewhere73. In addition, we included 401 individuals (including 25 neurotypical individuals) from the EU-AIMS LEAP74 and Paris75 cohorts. Across all cohorts, we included only unrelated individuals, who were predominantly of European Ancestry as defined by genetic principal components (5 SD deviations above or below the mean of PC1 and PC2 from the HapMap CEU population).

Additionally, we also included data from 1981 unrelated individuals (1000 males, 981 females) from the Nijmegen Biomedical Study (NBS) to provide support for the independent SNPs with P < 1 × 10−6 in the non-stratified GWAS. Participants were asked the question: ‘It upsets me if my daily routine is disturbed’, which is related to a non-social domain of autism, and is similar to an item in the Autism Spectrum Quotient. Further information including genotyping and quality control is provided elsewhere43. Genetic association for the top SNPs were conducted using age, sex, and the first five genetic principal components as covariates using linear regression.

Phenotypes

The primary phenotype for this study is the SQ-R, which was used to conduct a GWAS in participants from 23andMe. The SQ-R is self-report measure of systemising drive, or interest in rule-based patterns10. The SQ-R taps a variety of domains of systemising, such as interest in mechanical (e.g. car engines), abstract (e.g. mathematics), natural (e.g. the weather), motor (e.g. knitting), and collectible (e.g. stamp collecting) systems. There are 75 items on the SQ-R, with a maximum score of 150 and a minimum score of 0. Scores on the test are normally distributed10. The SQ-R has good cross-cultural stability and good psychometric properties with Cronbach’s alpha ranging from 0.79 to 0.94 in different studies76. Test–retest reliability available in a Dutch sample indicated high reliability of 0.79 (Pearson correlation)76. This was supported by another study in 4058 individuals which identified high internal cohesion77. Exploratory followed by confirmatory factor analysis using Rasch modelling suggests that the SQ-R is unidimensional77. A sex difference has been observed in multiple studies with males, on average, scoring significantly higher than females10,51. Criterion validity shows that the SQ-R has a modest but significant correlation with the mental rotation test (r = 0.25, P = 0.013), as well as its subscales78. Autistic individuals, on average, score higher on the SQ-R in multiple different studies10,51,79. Further, the SQ-R also predicts autistic traits, with a combination of the SQ-R and the Empathy Quotient predicting as much as 75% of the variance on the autism spectrum quotient, a measure of autistic traits10. The SQ-R has been validated using a short form in a very large population of 600,000 controls and 36,000 autistic individuals12.

In addition, we used the following secondary phenotypes: SCDC in ALSPAC, ADOS-G social and communication scores and the RBS-R in the other cohorts. We also used a single question which is a measure of ‘insistence on sameness’ in the NBS cohort.

The SCDC is a questionnaire that measures difficulties in verbal and nonverbal communication, and social interaction including reciprocal social interaction61. The questionnaire consists of 12 questions, with scores ranging from 0 to 24, and with higher scores reflecting difficulties in social interaction and communication. The SCDC has good internal consistency (0.93) and good test–retest reliability (0.81)61. The SCDC has reasonable specificity and sensitivity in distinguishing autistic from control individuals80. Previous research has demonstrated that the SCDC is genetically correlated with autism44,45,57. We conducted a GWAS of SCDC to investigate whether it is genetically correlated with SQ-R in this study. We used mother-reported SCDC scores on children aged 8. Although SCDC has been measured at different ages in the ALSPAC cohort, we chose SCDC scores measured at age 8 as this has the largest sample size and has high \({{h}}_{{\mathrm{SNP}}}^2\)19 (h2 = 0.24 ± 0.07).

We chose two measures of social and non-social traits in autistic individuals. For the social trait, we used the social and communication domain scores from the ADOS-G, a widely used instrument for diagnosing and assessing autism in four cohorts (SSC, AGRE, EU-AIMS LEAP, and Paris). Participants completed one of the following ADOS-G modules81: 1 (used for children with little or no phrase speech), 2 (for children with non-fluent speech), 3 (verbally fluent children), and 4 (verbally fluent adolescents and adults). For this study, we used the raw totals of the scores from the social domain and the communication domain, combined. Scores for all four modules range from 0 to 24. The ADOS-G has high overall internal consistency, and high test–retest reliability for the social and communication subscales81. The choice for combining the social and communication domain scores were informed by factor analysis which suggested that the two domains contribute to one underlying factor82.

In contrast to the social and communication domain, the restricted and repetitive behaviour domain of the ADOS-G has poor test–retest reliability (r < 0.6) and a smaller range of scores (0–8) as it captures fewer repetitive and restrictive behaviour81. Hence, for this study, we used sores on the RBS-R83. The RBS-R is a measure developed to specifically measure restricted and repetitive behaviours in autistic individuals and captures stereotyped, self-injurious, sameness, compulsive, ritualistic, and restricted behaviour84, and has high inter-rater reliability and internal consistency84. The RBS-R comprises 43 questions with scores ranging from 0 to 3 for each item based on a Likert scale.

‘Insistence on sameness’ in the NBS cohort was measured using a single item: ‘It upsets me if my daily routine is disturbed’. This is related to a non-social domain of autism, and is again similar to an item in the Autism Spectrum Quotient. Participants were asked to indicate on a 4-point Likert scale ‘definitely agree’, ‘slightly agree’, ‘slightly disagree’, ‘definitely disagree’.

Genotyping, imputation, and quality control and genetic association in the 23andMe cohort

Details of genotyping, imputation and quality control in the 23andMe cohort are provided elsewhere47. Briefly, unrelated participants were included if they had a call rate of >98.5%, and were of primarily European ancestry (97% European ancestry). A total of 1,030,430 SNPs (including InDels) were genotyped or imputed. SNPs were excluded if: they failed the Hardy–Weinberg equilibrium test at P < 10−20; had a genotype rate of <90%; they failed the parent-offspring transmission test using trio data in the larger 23andMe research participant database; or if allele frequencies were significantly different from the European 1000 Genomes reference data (χ2 test, P < 10−20). Phasing was conducted using Beagle (version 3.3.1)85 in batches of 8000–9000 individuals. This was followed by imputation against all-ethnicity 1000 Genomes haplotypes (excluding monomorphic and singleton sites) using Minimac286. Genetic association analyses were restricted to SNPs with a minor-allele frequency > 1%. After quality control, 9,955,952 SNPs (imputed and genotyped) were included in the GWAS.

Our primary analysis was an additive model of genetic effects and was conducted using a linear regression with age, sex, and the first five ancestry principal components included as covariates. In addition, given the modest sex difference, we also conducted sex-stratified analyses. SNPs were considered significant at a genome-wide threshold of P < 5 × 10−8. Leading SNPs were identified after LD-pruning using Plink (r2 > 0.8). Winner’s curse correction was conducted using an FDR-based shrinking87.

We calculated variance explained by first standardising the regression estimates and then squaring the estimates. This is equivalent to:

$$R^2 = \hat B_j^2\frac{{2\left( {{\mathrm{MAF}}_j} \right)(1 - {\mathrm{MAF}}_j)}}{{\sigma _y^2}},$$

where R2 is the proportion of variance explained for SNP j . \(\hat B_j^2\) is the non-standardised regression coefficient, MAF is the minor-allele frequency for SNP j , and \(\sigma _y^2\) is the variance of SQ. Further details of this formula are provided in the Supplementary Note.

Genotyping, imputation, and quality control and genetic association in the ALSPAC

The SCDC61 scores were calculated from children of the 90 s (ALSPAC cohort)71, in children aged 8. In total, SCDC scores were available on N = 7,825 children. From this, we removed individuals for whom complete SCDC scores were not available. After excluding related individuals and individuals with no genetic data, data was available on a total of N = 5,421 unrelated individuals.

Participants were genotyped using the Illumina® HumanHap550 quad chip by Sample Logistics and Genotyping Facilities at the Wellcome Sanger Institute and LabCorp (Laboratory Corportation of America) using support from 23andMe. Individuals were excluded based on gender mismatches, high missingness (>3%), and disproportionate heterozygosity. We restricted subsequent analyses to individuals of European descent (CEU), which were identified by multi-dimensional scaling analysis and compared with Hapmap II (release 22). Individuals were also removed if cryptic relatedness, assessed using identity by descent, was >0.1. Genotyped SNPs were filtered out if they had >5% missingness, violated Hardy–Weinberg equilibrium (P < 1 × 10−6), and had a minor-allele frequency < 1%, resulting in a total of 526,688 genotyped SNPs. Haplotypes were estimated using data from mothers and children using ShapeIT (v2.r644)88. Imputation was performed using Impute2 V2.2.289 against the 1000 genomes reference panel (Phase 1, Version 3). Imputed SNPs were excluded from all further analyses if they had a minor-allele frequency < 1% and info < 0.8. After quality control, there were 8,282,911 genotyped and imputed SNPs that were included in subsequent analyses.

Dosage data from BGEN files were converted using hard-calls, with calls with uncertainty > 0.1 treated as missing data. Post-imputation, we excluded SNPs that deviated from Hardy–Weinberg equilibrium (P < 1 × 10−6), with minor-allele frequency < 0.01 and missing call rates > 2%. We further excluded individuals with genotype missing rates > 5%. The SCDC score was not normally distributed so we log-transformed the scores and ran regression analyses using the first two ancestry principal components and sex as the covariates using Plink 2.0 (ref. 90).

The log-transformed SCDC scores (henceforth, SCDC scores) had a modest but significant \({{h}}_{{\mathrm{SNP}}}^2\) as quantified using LDSR \(\left( {{{h}}_{{\mathrm{SNP}}}^2 = 0.12 \pm 0.05} \right)\). LDSR intercept (0.99) suggested that there was no inflation in GWAS estimates due to population stratification. The λ GC was 1.013. We replicated the previously identified genetic correlation with autism57 (constrained intercept) using our SCDC GWAS (r g = 0.45 ± 0.18, P = 0.01). In addition, we also identified a negative genetic correlation between educational attainment58 and SCDC (r g = −0.30 ± 0.11, P = 0.007).

Genomic inflation factor, heritability, and functional enrichment for the SQ-R GWAS

LDSR91,92 was used to calculate for inflation in test statistics due to unaccounted population stratification. Heritability was calculated using LDSR using the north-west European LD scores. Difference in heritability between males and females was quantified using:

$$Z = \frac{{h_{{\mathrm{males}}}^2 - h_{{\mathrm{females}}}^2}}{{\sqrt {{\mathrm{SE}}_{{\mathrm{males}}}^2 + {\mathrm{SE}}_{{\mathrm{females}}}^2} }},$$

where Z is the Z-score for the difference in heritability for a trait, \(\left( {{{h}}_{{\mathrm{males}}}^2 - {{h}}_{{\mathrm{females}}}^2} \right)\) is the difference \({{h}}_{{\mathrm{SNP}}}^2\) estimate in males and females, and SE is the standard errors for heritability. Two-tailed P-values were calculated and reported as significant if P < 0.05.

For the primary GWAS (non-stratified analyses), we conducted functional annotation using FUMA93. We restricted our analyses to the non-stratified analyses due to the high genetic correlation between the sexes and the low statistical power of the sex-stratified GWAS. We conducted gene-based association analyses using MAGMA94 within FUMA and report significant genes after using a stringent Bonferroni corrected P < 0.05. In addition, we conducted enrichment for tissue specific expression and pathway analyses within FUMA. For the significant SNPs, we investigated enrichment for eQTLs using brain tissues in the BRAINEAC and GTEx95 database within FUMA. We further conducted partitioned heritability for tissue-specific active chromatin marks and baseline functional categories using extended methods in LDSR96.

Hi-C-based annotations of fine mapped loci

We fine mapped three genome-wide significant loci (index SNPs: rs4146336 and rs1559586 or SQ; rs8005092 for SQ-R males) to obtain credible SNPs. First, we selected SNPs with P < 0.01 that are located in the LD region (r2 > 0.6) with an index SNP. LD structure within a locus was constructed by calculating correlations between SNPs within a locus (1KG v20130502). CAVIAR97 was then applied to the summary association statistics and LD structure for each index SNP to generate potentially causal (credible) SNPs with a posterior probability of 0.95. In total, we identified 14 credible SNPs from the three GWS loci.

For each locus, candidate genes were identified by mapping credible SNPs based on physical interactions in foetal brain as previously described98. One locus (index SNP rs4146336) was mapped to two genes, LSAMP and PTMAP8, indicating that two credible SNPs (rs13066948 and rs11713893) located in this locus physically interact with these genes.

Genetic correlation

For all phenotypes, we performed genetic correlation without constraining the intercept using LDSR. We identified significant genetic correlations using a Bonferroni adjusted P-value < 0.05. For the primary genetic correlation analysis with SQ-R, we included psychiatric conditions57,63,99,100,101,102, personality traits103,104,105, measures of intelligence58,59,106,107, and social traits related to autism46,55 including scores on the SCDC, as previous research has investigated the phenotypic correlation between these domains and systemising10,78,108,109,110,111,112.

To understand the correlation between systemising and various phenotypes that have been genetically correlated with autism, we used GWAS data from 15 phenotypes including autism. 10 of these phenotypes (cognitive aptitude59, educational attainment58, tiredness113, neuroticism103, subjective wellbeing103, schizophrenia114, major depression102, depressive symptoms103, ADHD63, and chronotype115), have been previously reported to be significantly genetically correlated with autism out of 234 phenotypes tested using LDHub62 (P < 2.1 × 10−4). We excluded college degree from this list, as previous work has identified near perfect genetic correlation between educational attainment and college degree58. In addition, we included data from friendship satisfaction55, family relationship satisfaction55, systemising, and self-reported empathy46, all of which are also significantly genetically correlated with autism with P < 2.1 × 10−4. These four additional phenotypes were not included in the previous paper which investigated genetic correlations with autism. Details of sample sizes with PMIDs/DOIs are provided in Supplementary Data 12. Cross trait genetic correlations were computed for all 15 phenotypes, and results were corrected for multiple testing using Bonferroni correction. A correlogram was created after using hierarchical clustering to cluster the phenotypes.

To investigate whether the combination of negative genetic correlation social traits and positive genetic correlation for non-social traits is specific to autism, we conducted a genetic correlation between all psychiatric conditions for which we had access to summary GWAS statistics (ADHD63, Anxiety116, Autism57, Anorexia101, Bipolar Disorder99, Major Depressive Disorder102, OCD117,118, PTSD119, and Schizophrenia114) and SQ-R, self-reported empathy measured using the EQ46 and friendship satisfaction55. We chose friendship satisfaction and self-reported empathy as representative of social traits as these are the most relevant to the social domain of autism for which we had access to GWAS summary statistics. The EQ is a short, 40-item self-report measure of empathy, which has been widely used and has good psychometric properties60,120. In addition, differences in aspects of empathy compared to the neurotypical population have been widely reported in autism50,51,121, and empathy is one of the items in measures such as ADOS-G. Friendship satisfaction was selected as difficulties in making friends is listed as a criteria for an autism diagnosis in the DSM-51.

GWIS and GSEM

To investigate whether the SQ-R is genetically correlated with autism independent of the genetic effects of educational attainment, we constructed a unique SQ-R phenotype after conditioning on the genetic effects of educational attainment using GWIS122. GWIS takes into account the genetic covariance between the two phenotypes to calculate the unique component of the phenotypes as a function of the genetic covariance and the \({{h}}_{{\mathrm{SNP}}}^2\). Before performing GWIS, we standardised the beta coefficients for the SQ-R GWAS by using the following formula:

$$\widehat {B_{{\mathrm{std}}}} = \hat B\sqrt {\frac{{2\left( {{\mathrm{MAF}}} \right)\left( {1 - {\mathrm{MAF}}} \right)}}{{\sigma _y^2}}},$$

where \(\widehat {B_{{\mathrm{std}}}}\) is the standardised regression coefficients, \(\hat B\) is the regression coefficient obtained from the non-standardised GWAS, MAF is the minor-allele frequency, \(\sigma _y^2\) is the variance of the SQ-R. This equation is explained in detail in the Supplementary Note. We conducted GWIS using only educational attainment as we were unclear whether the GWAS of cognitive aptitude59 was conducted on a standardised phenotype. Further, there is a high genetic correlation between cognitive aptitude and educational attainment. In addition to GWIS, to validate the findings, we conducted GSEM123, a complementary but independent method. GSEM uses the genetic correlations and covariances calculated using LDSR after accounting for sample overlap.

Polygenic scores in the SSC, AGRE, EU-AIMS LEAP, and Paris cohorts

We generated polygenic scores for SQ-R (mean weighted score of all the alleles that contribute to higher systemising) in 2221 probands from the SSC (Discovery dataset). We downloaded genotype data from the SSC from SFARI base (https://www.sfari.org/resource/sfari-base/). Individuals were genotyped on three different platforms: Illumina Omni2.5, Illumina 1Mv3, or Illumina 1Mv1. Informed consent or assent was obtained from all participants. In addition, the research team obtained ethical approval from the Cambridge Human Biology Research Ethics Committee to access and analyse the de-identified data from the SSC. We conducted stringent quality control and imputation separately for each platform. The full pipeline is available here: https://github.com/autism-research-centre/SSC_liftover_imputation. Briefly, individuals were excluded if they had: a genotyping rate < 95%, excessive or low heterozygosity (less or more than 3 SD from the mean), mismatched reported and genetic sex, and families with Mendelian errors > 5%. We further removed SNPs that significantly deviated from Hardy–Weinberg equilibrium (P < 1 × 10−6), had Mendelian errors in >10% of the families, and SNPs that were not genotyped in >10% of the families. We then conducted multi-dimensional scaling using the HapMap3 phase 3 population using the unrelated individuals CEU and TSI populations as representatives of the European population. This was conducted only in the parents to retain unrelated individuals for multi-dimensional scaling. Genetic principal components were calculated using only SNPs with minor-allele frequency > 5%, and pruning the SNPs in Plink using an r2 of 0.2. We excluded families from further downstream analyses if either one the parents were greater or less than 5 standard deviations from the means of the first two genetic principal components calculated using only the unrelated individuals in HapMap3 CEU and TSI populations. Quality control was done using Plink v 1.9 and R. Phasing and imputation were conducted using the Michigan Imputation Server (https://imputationserver.sph.umich.edu/start.html) using the 1000 genomes Phase 3 v5 as the reference panel.

Polygenic scores were generated using PRSice2 (https://choishingwan.github.io/PRSice/) for the SQ-R using the non-stratified GWAS data. We calculated the mean polygenic score for each of the 2221 probands in the SSC, after clumping SNPs using an R2 threshold of 0.1. Prior to generating polygenic scores, we confirmed that the probands were not related to each other using identity by descent PI-HAT > 0.15 as a relatedness cut-off. We used a P-value threshold of 1 as previous research on educational attainment, subjective wellbeing and social relationship satisfaction, all suggest that the maximum variance explained is at a threshold of 1 (refs. 58,103). This is expected for highly polygenic traits where many SNPs incrementally contribute to the variance explained124. Polygenic scoring was done using standardised scores on two different phenotypes as the dependent variable (RBS-R and the social and communication domain of the ADOS-G). We included sex, platform, the first 15 genetic principal components and standardised full-scale IQ as covariates. In addition, for the analysis of ADOS-G, we included the ADOS-G module as a covariate. Linear regression was conducted in R. A total of 135,233 SNPs were included in the polygenic score analyses after clumping and thresholding.

To validate the polygenic scores, we conducted additional polygenic score analysis using data combined from the AGRE, EU-AIMS LEAP and Paris cohorts. We followed similar quality control and imputation procedures to the SSC cohort. Given that this dataset was a mix of related and unrelated individuals, we chose unrelated individuals using a genomic relationship matrix (GRM) as provided in GCTA (grm-cutoff 0.05)125. To calculate GRMs, we included only SNPs with minor-allele frequency > 1%. Scripts are provided here: https://github.com/vwarrier/PARIS_LEAP_analysis. Polygenic scores were calculated using PRSice2 as described for the SSC data. Given the differences in dataset, polygenic scores were calculated separately for the AGRE dataset, and the EU-AIMS LEAP and Paris datasets combined. For each regression, we included sex and the first ten genetic principal components (standardised). The dependent variables were standardised scores on the RBS-R (N = 426) and the ADOS-G social and communication subscale (N = 475). IQ information was unavailable for most individuals, and hence we did not include IQ as a covariate. We combined the results of the EU-AIMS LEAP and Paris cohorts, and the AGRE dataset using inverse-variance weighted fixed-effect meta-analysis using the formula below:

$$\begin{array}{*{20}{l}} {w_i} \hfill & = \hfill & {1/{\mathrm{SE}}_i^2} \hfill \\ {{\mathrm{SE}}_{{\mathrm{meta}}}} \hfill & = \hfill & {\sqrt {1/{\mathrm{\Sigma }}_iw_i} } \hfill \\ {Beta_{{\mathrm{meta}}}} \hfill & = \hfill & {{\mathrm{\Sigma }}_i\beta _iw_i/{\mathrm{\Sigma }}_iw_i} \hfill \end{array},$$

where β i is the standardised regression coefficient of the polygenic scores, SE i is the associated standard error, and w i is the weight.

Bivariate GREML

We conducted bivariate genetic correlation using GCTA GREML to test the genetic correlation between the ADOS social and communication domains and the RBS-R scores. We created a GRM after including autistic individuals from the SSC, AGRE, EU-AIMS LEAP, and Paris cohorts. We excluded SNPs and individuals using the same quality control pipeline as applied to the SSC dataset outlined in the section above. We further restricted our analysis only to SNPs with a minor-allele frequency > 1%. We excluded related individuals (–grm-cutoff 0.05) resulting in a total of 2989 individuals. Of this, 2652 individuals had scores for the ADOS social and communication domain and 2550 individuals had scores on the RBS-R. We included sex and the first ten genetic principal components as covariates.

Reporting summary

Further information on research design is available in the Nature Research Reporting Summary linked to this article.