A 3D approach to understanding heterogeneity in early developing autisms

Abstract

Background: Phenotypic heterogeneity in early language, intellectual, motor, and adaptive functioning (LIMA) features are amongst the most striking features that distinguish different types of autistic individuals. Yet the current diagnostic criteria uses a single label of autism and implicitly emphasizes what individuals have in common as core social-communicative and restricted repetitive behavior difficulties. Subtype labels based on the non-core LIMA features may help to more meaningfully distinguish types of autisms with differing developmental paths and differential underlying biology.

Methods: Unsupervised data-driven subtypes were identified using stability-based relative clustering validation on publicly available Mullen Scales of Early Learning (MSEL) and Vineland Adaptive Behavior Scales (VABS) data (n = 615; age = 24-68 months) from the National Institute of Mental Health Data Archive (NDA). Differential developmental trajectories between subtypes were tested on longitudinal data from NDA and from an independent in-house dataset from UCSD. A subset of the UCSD dataset was also tested for subtype differences in functional and structural neuroimaging phenotypes and relationships with blood gene expression. The current subtyping model was also compared to early language outcome subtypes derived from past work.

Results: Two autism subtypes can be identified based on early phenotypic LIMA features. These data-driven subtypes are robust in the population and can be identified in independent data with 98% accuracy. The subtypes can be described as Type I versus Type II autisms differentiated by relatively high versus low scores on LIMA features. These two types of autisms are also distinguished by different developmental trajectories over the first decade of life. Finally, these two types of autisms reveal striking differences in functional and structural neuroimaging phenotypes and their relationships with gene expression and may highlight unique biological mechanisms.

Limitations: Sample sizes for the neuroimaging and gene expression dataset are relatively small and require further independent replication. The current work is also limited to subtyping based on MSEL and VABS phenotypic measures.

Conclusions: This work emphasizes the potential importance of stratifying autism by a Type I versus Type II distinction focused on LIMA features and which may be of high prognostic and biological significance.

Keywords: Clustering; Gene expression; Stratification; Subtypes; fMRI.

Fig. 1

Schematic of data analysis workflow/pipeline. This figure shows a workflow of the analyses conducted in this work. Panels A and B depict how we select, clean, and preprocess data from the National Institute of Mental Health Data Archive (NDA). To identify subtypes in a data-driven manner, we used stability-based relative clustering validation analysis (reval) (panel C) on a final sample of 615 autistic children aged 24–68 months. Subtypes are tested for their prognostic validity (panel D) utilizing longitudinal data from NDA. The stratification model was tested again for developmental trajectory differences in an independent longitudinal dataset from UCSD ACE (panels E–F). Subtypes were also tested for neurobiological differentiation using structural and functional neuroimaging data and relationships to blood leukocyte gene expression patterns (panel G)

Fig. 2

Robust and highly generalizable early distinctions between Type I versus Type II autisms. Panel A shows that a 2-subtype solution is the optimal clustering solution that minimizes normalized stability via stability-based relative clustering validation (reval). This solution generalizes with 98% accuracy in the held-out NDA validation dataset. Panel B shows UMAP density plots that depict separated peaks of distributions for Type I (blue) and Type II (orange) autism subtypes. Panel C plots the standardized scores across all MSEL and VABS features to describe how subtypes manifest as relatively low (Type I; blue) versus high (Type II; orange) scores across all features. Panels D-E show plots of developmental trajectories for the Type I (blue) and Type II (orange) autisms, respectively for NDA (D) and UCSD ACE (E) datasets

Fig. 3

Functional neural response to speech and associations with gene expression. Panel A depicts whole-brain activation analysis maps for each of the groups under study. All maps are thresholded at voxel-FDR q < 0.05, except for LD/DD, Poor and Type I groups, which are thresholded at cluster-forming p < 0.05 and cluster FDR q < 0.05 for visualization purposes. Panel B depicts standardized effect sizes (Cohen’s d) for the LH temporal ROI for all pairwise group comparisons for A3D (top) and ELO (bottom) stratification models. The asterisks and green outlined cells indicate statistical significance at p < 0.05. Similar effect size heatmaps are shown in panel C for the RH temporal ROI. Panels D and E show results for the PLS analysis of the statistically significant LV1 fMRI-gene expression relationship identified in the ELO (D) and A3D (E) models. The whole-brain map in these panels show brain bootstrap ratios (BSR) that indicate which voxels most strongly contribute to the fMRI-gene expression LV1 relationship. The heatmaps to the right on each panel shows each gene co-expression module along the rows and groups along the columns. The modules of importance are highlighted with green outlines (non-zero modules) as they contain a correlation with fMRI response to speech whose 95% bootstrap CIs do not encompass a correlation of zero

Fig. 4

Gene expression associations with cortical surface area. Results for PLS analyses on the phenotype of cortical surface area for the ELO (A) or A3D (B, A3D SA LV2; C, A3D SA LV1) models. The parcellated brain maps show brain bootstrap ratios (BSR), whereby more extreme BSRs indicate brain regions that more strongly contribute to the LV relationship. The heatmap on the right of each panel shows correlations with each gene co-expression module along the rows and each group along the columns. The important non-zero modules are highlighted as cells with green outlines

Fig. 5

Gene expression associations with cortical thickness. Results for PLS analyses on the phenotype of cortical thickness for the ELO (A, ELO CT LV1; C, ELO CT LV2) or A3D (B, A3D CT LV1) models. The parcellated brain maps show brain bootstrap ratios (BSR), whereby more extreme BSRs indicate brain regions that more strongly contribute to the LV relationship. The heatmap on the right of each panel shows correlations with each gene co-expression module along the rows and each group along the columns. The important non-zero modules are highlighted as cells with green outlines

Fig. 6

Enrichment analysis showing overlap between attenuated transcriptomic regional identity (ARI) genes detected in post-mortem cortical tissue and non-zero modules detected in imaging-gene expression PLS analyses. Panel A depicts a schematic of how these enrichment analyses were conducted. A gene set known as attenuated transcriptomic regional identity (ARI) genes were extracted from a prior study on post-mortem cortical tissue in autistic patients [42]. These ARI genes are genes that have substantial between-region gene expression differences in non-autistic brains, but much more attenuated regional identity differences in gene expression in autism. ARI genes are depicted in blue in the Venn diagram. In green, the Venn diagram shows our set of genes isolated from PLS imaging-gene expression association analyses (i.e. genes from non-zero modules). The degree of overlap is then tested with enrichment odds ratios and hypergeometric p values from gene set enrichment analyses. Panel B shows a heatmap of the enrichment odds ratios (numbers in the center of each cell) and where color in each cell indicates the -log10(p-value) for each hypergeometric test. Cells with a thick black outline indicate enrichment tests that pass at FDR q < 0.05, while cells with smaller thin black outline pass at FDR q < 0.1. The columns in the heatmap represent comparisons when the gene list derives from the ELO (left) or A3D (right) model. Each row of the heatmap indicates a different gene list extracted from PLS analyses, with the top row indicating genes from non-zero modules in PLS analyses that examine fMRI responses to speech. The remaining rows depict comparisons where the PLS gene lists come from associations with surface area (SA) or cortical thickness (CT) and we have annotated which rows can be interpreted as effects indicative of ‘normative patterning’ effects (effects driven by TD correlations and which are shared in Good ELO or Type II A3D subtypes) versus atypical patterning effects (e.g., effects driven specifically by the ELO Poor or Type I A3D subtypes)