Structural models of genome-wide covariance identify multiple common dimensions in autism

Abstract

Common genetic variation has been associated with multiple phenotypic features in Autism Spectrum Disorder (ASD). However, our knowledge of shared genetic factor structures contributing to this highly heterogeneous phenotypic spectrum is limited. Here, we developed and implemented a structural equation modelling framework to directly model genomic covariance across core and non-core ASD phenotypes, studying autistic individuals of European descent with a case-only design. We identified three independent genetic factors most strongly linked to language performance, behaviour and developmental motor delay, respectively, studying an autism community sample (N = 5331). The three-factorial structure was largely confirmed in independent ASD-simplex families (N = 1946), although we uncovered, in addition, simplex-specific genetic overlap between behaviour and language phenotypes. Multivariate models across cohorts revealed novel associations, including links between language and early mastering of self-feeding. Thus, the common genetic architecture in ASD is multi-dimensional with overarching genetic factors contributing, in combination with ascertainment-specific patterns, to phenotypic heterogeneity.

Fig. 1. Workflow of the study.

a Multi-stage study design. Multivariate discovery analyses were carried out in the Simons Powering Autism Research (SPARK) sample (Stages I-III) and the best-fitting model in SPARK was followed-up in the Simons Simplex Collection (SSC, Stage IV). b Data-driven genomic covariance modelling approach, including a step-wise combination of principal component analysis (PCA), exploratory factor analysis (EFA) and Genetic-relationship-matrix structural equation modelling (GRM-SEM), as described in the Methods.

Fig. 2. GREML heritability estimates for SPARK and SSC phenotypes.

a GREML h²_SNP of continuous and categorical ASD phenotypes with p ≤ 0.1 in the SPARK sample (N ≤ 5132). A complete figure of all analysed phenotypes is shown in Supplementary Fig. 2. Information on phenotype description, sample size and exact heritability and p-values is available in Supplementary Data 1. b GREML h²_SNP of continuous and categorical ASD phenotypes in the SSC sample (N ≤ 1940). Information on phenotype description, sample size and exact heritability and p-values is available in Supplementary Data 6. The error bars represent standard errors. Evidence for GREML h²_SNP estimates was based on likelihood ratio tests. No adjustments for multiple-testing were carried out. Estimates were based on transformed scores: deviance residuals (for categorical phenotypes) or rank-transformed residuals (for continuous phenotypes). DCDQ (Developmental Coordination Disorder Questionnaire), GREML (Genome-based restricted maximum likelihood), h²_SNP (Single nucleotide polymorphism-based heritability), ODD (oppositional defiant disorder), RBSR (Repetitive Behaviour Scale-Revised).

Fig. 3. Best-fitting model in SPARK.

a Scree plot based on the eigenvalue decomposition of genetic correlations derived from a GRM-SEM Cholesky model, depicting the number of estimated shared genetic factors (in black) according to an optimal coordinate criterion. The dashed line indicates the “scree” of the plot (grey). b Path diagram depicting the best-fitting multi-factor GRM-SEM IPC model. Observed measures are represented by squares and latent variables by circles (A: shared genetic factor, AS: specific genetic factor, E: residual factor). Single-headed arrows define factor loadings (shown with their corresponding SEs). The genetic part of the model has been modelled using an Independent Pathway model. Grey dotted and coloured solid arrows define shared genetic factor loadings with p > 0.05 and p ≤ 0.05, respectively. Black dotted lines define specific genetic factor loadings with p > 0.05. The residual part has been modelled using a Cholesky model and all residual factor loadings are shown in grey. The full parameter table is shown in Supplementary Data 5. Evidence for GRM-SEM factor loadings was assessed with Wald tests (two-sided). Given the multivariate design, no adjustment for multiple comparisons was carried out. c Corresponding standardised genetic variance (GRM-SEM h²_SNP) plot. SEs for GRM-SEM h²_SNP contributions have been omitted for clarity. d Corresponding correlogram of genetic correlations (r_g). Numeric values for genetic correlations that are not predicted by the genetic model structure were omitted. A_lang (genetic language performance factor), A_dev (genetic developmental motor delay factor), A_beh (genetic behavioural-problems factor), DCDQ (Developmental Coordination Disorder Questionnaire), h²_SNP (Single nucleotide polymorphism-based heritability), IPC (Independent Pathway-Cholesky GRM-SEM model), ODD (Oppositional Defiant Disorder), RBSR (Repetitive Behaviours Scale-Revised).

Fig. 4. Follow-up multi-factor GRM-SEM model in the SSC.

a Scree plot based on the eigenvalue decomposition of genetic correlations derived from a GRM-SEM Cholesky model, depicting the number of estimated shared genetic factors (in black) according to an optimal coordinate criterion. The dashed line indicates the “scree” of the plot (grey). b Path diagram depicting the best-fitting multi-factor GRM-SEM IPC model based on largely comparable phenotypes as studied in SPARK. Observed measures are represented by squares and latent variables by circles (AF: shared genetic factor, AS: specific genetic factor, E: residual factor). Single-headed arrows define factor loadings (shown with their corresponding SEs). The genetic part of the model has been modelled using an Independent Pathway model. Grey dotted and coloured solid arrows define shared genetic factor loadings with p > 0.05 and p ≤ 0.05, respectively. Black dotted lines define specific genetic factor loadings with p > 0.05. The residual part has been modelled using a Cholesky model and all residual factor loadings are shown in grey. The full parameter table is shown in Supplementary Data 7. Evidence for GRM-SEM factor loadings was assessed with Wald tests (two-sided). Given the multivariate design, no adjustment for multiple comparisons was carried out. c Corresponding standardised genetic variance (GRM-SEM h²_SNP) plot. SEs for GRM-SEM h²_SNP contributions have been omitted for clarity. d Corresponding correlogram of genetic correlations (r_g). Numeric values for genetic correlations that are not predicted by the genetic model structure were omitted. A_F1,2,3 (Genetic factor 1,2,3), h²_SNP (Single nucleotide polymorphism-based heritability), IPC (Independent Pathway-Cholesky GRM-SEM model), ODD (Oppositional Defiant Disorder).

Fig. 5. Characterisation of identified genetic factor structures in SPARK.

a Path diagram of an extended GRM-SEM IPC model mapping liability to Asperger (reference: Asperger against other ASD subcategories) onto the model structure of the best-fitting model in SPARK. b Corresponding standardised genetic variance (GRM-SEM h²_SNP) plot. SEs for GRM-SEM h²_SNP contributions have been omitted for clarity. c Genetic correlations with liability to Asperger. d Path diagram of an extended GRM-SEM IPC model mapping polygenic scores for educational attainment (PGS_EA) onto the model structure of the best-fitting model in SPARK. e Corresponding standardised genetic variance (GRM-SEM h²_SNP) plot. SEs for GRM-SEM h²_SNP contributions have been omitted for clarity. f Genetic correlations with PGS_EA. a, d Observed measures are represented by squares and latent variables by circles (A_lang/A_dev/A_beh: shared genetic factor, AS: specific genetic factor, E: residual factor). Single-headed arrows define factor loadings (shown with their corresponding SEs). The genetic part of the model has been modelled using an Independent Pathway model. Grey dotted and coloured solid arrows define shared genetic factor loadings with p > 0.05 and p ≤ 0.05, respectively. Black dotted lines define specific genetic factor loadings with p > 0.05. Factor loadings for the mapping variable are shown in blue (dotted: p > 0.05; solid p ≤ 0.05). The residual part has been modelled using a Cholesky model (grey). Evidence for GRM-SEM factor loadings was assessed with Wald tests (two-sided). Given the multivariate design, no adjustment for multiple comparisons was carried out. A_lang (genetic language performance factor), A_dev (genetic developmental motor delay factor), A_beh (genetic behavioural-problems factor), DCDQ (Developmental Coordination Disorder Questionnaire), h²_SNP (single nucleotide polymorphism-based heritability), IPC (Independent Pathway-Cholesky GRM-SEM model), ODD (Oppositional Defiant Disorder), RBSR (Repetitive Behaviours Scale-Revised), r_g (genetic correlation).