Gray matter covariations and core symptoms of autism: the EU-AIMS Longitudinal European Autism Project

Background Voxel-based morphometry (VBM) studies in autism spectrum disorder (autism) have yielded diverging results. This might partly be attributed to structural alterations being associating with the combined influence of several regions rather than with a single region. Further, these structural covariation differences may relate to continuous measures of autism rather than with categorical case–control contrasts. The current study aimed to identify structural covariation alterations in autism, and assessed canonical correlations between brain covariation patterns and core autism symptoms. Methods We studied 347 individuals with autism and 252 typically developing individuals, aged between 6 and 30 years, who have been deeply phenotyped in the Longitudinal European Autism Project. All participants’ VBM maps were decomposed into spatially independent components using independent component analysis. A generalized linear model (GLM) was used to examine case–control differences. Next, canonical correlation analysis (CCA) was performed to separately explore the integrated effects between all the brain sources of gray matter variation and two sets of core autism symptoms. Results GLM analyses showed significant case–control differences for two independent components. The first component was primarily associated with decreased density of bilateral insula, inferior frontal gyrus, orbitofrontal cortex, and increased density of caudate nucleus in the autism group relative to typically developing individuals. The second component was related to decreased densities of the bilateral amygdala, hippocampus, and parahippocampal gyrus in the autism group relative to typically developing individuals. The CCA results showed significant correlations between components that involved variation of thalamus, putamen, precentral gyrus, frontal, parietal, and occipital lobes, and the cerebellum, and repetitive, rigid and stereotyped behaviors and abnormal sensory behaviors in autism individuals. Limitations Only 55.9% of the participants with autism had complete questionnaire data on continuous parent-reported symptom measures. Conclusions Covaried areas associated with autism diagnosis and/or symptoms are scattered across the whole brain and include the limbic system, basal ganglia, thalamus, cerebellum, precentral gyrus, and parts of the frontal, parietal, and occipital lobes. Some of these areas potentially subserve social-communicative behavior, whereas others may underpin sensory processing and integration, and motor behavior.


FSIQ
To validate the results not being biased by low IQ participants we excluded the participants in Schedule D (FSIQ < 75) and performed an analogous 100-dimensional ICA decomposition, and automatic dimensional estimated factorizations (85 ICs) followed by post-hoc statistics again. In the 100-dimensional factorization case, we found the components with significant group effect were equivalent to IC10 and IC14 in the original analysis with the Schedule D participants included ( Figure S1a), however, they did not survive FDR correction (p<1.600×10 -5 ). Similarly, in the automatic dimensional factorization case, we found significant ICs corresponding to IC10 and IC14, however the IC corresponding to IC14 did not survive FDR correction (p<7.531×10 -4 ) ( Figure S1b). Figure S1. The components showed significant case-control differences. Panel a shows the components in 100 dimensional factorization excluding Schedule D participants. IC9 corresponds to IC14 in the 100 dimensional factorization including Schedule D participants (p=0.001), and IC12 corresponds to IC10 (p=1.600×10 -5 ). Neither of these ICs survive FDR correction (p<1.600×10 -5 ). Panel b shows the components in automatic dimensional factorization excluding Schedule D participants. IC10 corresponds to IC14 in 100 dimensional factorization including Schedule D participants (p=7.531×10 -4 ), and IC13 (p=3.830×10 -5 ) corresponds to IC10 that survived FDR correction (p<7.531×10 -4 ). The component maps were thresholded at 3<|Z|<5. IC, independent component.

ADHD symptoms
We ran additional group effect analyses (autism vs control) with ADHD symptoms as a covariate by using a dimensional score of ADHD symptoms (instead of the dummy code). The dimensional score was assessed by subscales of ADHD rating scale for symptoms of inattention and hyperactivity/impulsivity (Table  S3). The result slightly differed from the analysis using the ADHD categorical score, where IC14 was found with significant group effect but IC10 was not. In the current analysis with ADHD dimensional score, component 14 was found with no significant group effect after FDR either (p=0.004). Since the continuous scores provide more variance of ADHD symptoms, this finding suggests IC14 (amygdala, hippocampus, and parahippocampal gyrus) is possibly confounded by the variance of comorbid ADHD symptoms.
Moreover, regarding the neural complexity of comorbidity with autism and ADHD, we excluded the participants that fulfilled the ADHD diagnosis on the ADHD rating scale and the participants with ADHD rating score unavailable, and then re-ran ICA and statistical analyses. The covariates in group effect analyses were the same as in the original main analyses. This resulted in 160 (46.1% of 347 participants) individuals with autism and 180 (71.4% of 252 participants) individuals with TD being included. Not surprisingly given the significant reduction in numbers, we didn't find any significant result in these analyses. We attribute the non-significant output probably to the reduced statistical power.

Age
As the rate of gray matter development may not be linear across a wide age range (1), potential effects of age squared, age-by-group, and age squared-by-group interactions should be accounted for in the statistical model. Accordingly, to avoid overfitting the model and acquire an optimal mode of the statistical model, we additionally ran a full model with these variables on the components found significant with group effect, and stepwise removed the age relating terms depending on their contribution to the model fit. After comparing the goodness of fits between the models, we hence acquired an appropriate model for case-control difference analyses. In the main results, we found that autism group showed significant difference of IC10 and IC14 from TD group. Therefore, we did the above mentioned analyses on these two components separately. We observed that in the analysis of IC10 the age squared variable did improve the fit of statistical model, while none of the age relating variables enhanced the fit of model for IC14 (Table S4). Consequently, we additionally added the age squared variable into the original model, and then found similar outputs as the main results. That is, we found IC10 (b=-0.147, p=8.996X10-5) and IC14 (b=-0.132, p=5.465x10-4) remaining significance contributors to the group effects (FDR corrected, p<7.751x10-4).

Sex-by-group interaction
As sex ratio is uneven between autism and control group, we accounted for the possible effect of sexby-group interaction. Therefore, we repeated the analyses similarly to the Age section to find an optimal mode of the statistical model. As a result, adding sex-by-group did not significantly improve the fit of the model (Table S5). Therefore, we kept the original model.

Anxiety and depression symptoms
In addition to ADHD rating scale, we used the Development and Well-Being Assessment (DAWBA) anxiety and depression prediction scores to investigate their separate effects on group differences of structural covariance (2). In DAWBA, each scale reflects six levels of predication (i.e., from ~0.1% to >70%) of the probability of meeting clinically relevant diagnostic criteria for a disorder. The anxiety prediction score reflects the highest risk of an individual across a group of anxiety disorders (obsessive-compulsive disorder (OCD), generalized anxiety, panic disorder, agoraphobia, PTSD, separation anxiety, social phobia, and specific phobia). The depression prediction score was generated for major depression. The information for the participants with available score was shown in Table S6.
We used anxiety and depression scores as additional covariates in the statistical model separately. In the analyses including the anxiety score (N=494), no significant group effects found, while the component with smallest p value is component 14 (p=7.289x10 -4 ), which was found with significant group effect in original analyses, comprising amygdala, hippocampus, and parahippocampal gyrus. In the other analyses including depression score (N=446), component 10 (insula, frontal areas, and caudate) was found significantly different in autism group (p=1.302x10 -4 , FDR corrected, p<8.99x10 -4 ).
Although there were no significant results remaining in the analysis where anxiety score was accounted for, the p values are relatively small, and it's meaningful that the two components demonstrate differences when taking the anxiety and depression comorbidities into account. Previous studies showed individuals with autism and individuals with anxiety both involved in the structural differences in inferior frontal gyrus, insula, striatum and amygdala (3,4), which indicates that IC10 probably reflects shared variances related to autism and anxiety on structural covariance. Major depression was formerly reported related to gray matter alterations in amygdala, hippocampus (e.g. (5)), which may suggest structural covariance in IC14 reflect shared variances between autism and depression symptoms.

Medication use
Since there is no strong and specific medication for autism, along with the high rate of comorbidity, medication use is heterogeneous in its nature (e.g. antidepressant, antipsychotic, sedatives, and medication treating ADHD). Therefore, we used a categorical score as a covariate to indicate whether the individuals take psychotropic medication (acting on the nervous system [data available on N=597, N=148 on psychotropic medication, Table S6). After regressing out medication use, we found that IC14 was not significant after FDR correction (p=0.003), while IC10 still showed significant group differences (p=1.162x10 -5 , FDR corrected, p<4.406x10 -4 ). This is partly in line with findings that suggest the volume of subcortical area is associated with medication use (6, 7).
Unfortunately, as unknown medication use could be confounding, and the various medications that individuals use are complex, the results of adding detailed medication use as a covariate would be difficult to interpret at best and likely underpowered to enable specific conclusions.

Sample homogeneity
Considering potential diagnostic difference of data quality that might influence the results, we checked the group differences of mean correlation from sample homogeneity measure while regressing out additional confounders; sex, site, FSIQ and age. Moreover, we also added it as an additional covariate into the statistical model for detecting group effect on structural covariance. We found that the homogeneity of gray matter images in autism group (mean: 0.878, standard deviation: 0.007) had no significant difference from TD group (mean: 0.879, standard deviation: 0.006; b=-0.051, p=0.176). Furthermore, as a covariate, the image quality had no significant effect on the main results we found. That is, IC10 was still found significantly associated with autism group (p=1.788x10 -4 ), while the group effect of IC14 was found at FDR threshold (p=0.001).

Reproducibility of CCA results
To assess the reproducibility of CCA results, we employed a leave-one-out (LOO) approach by randomly resampling subsets of the sample with participant number from 50 to 325 (ADI&ADOS)/194 (SRS&RBS&SSP), and repeating each LOO analysis 50 times. In each subset, we separately correlated the main mode weights of brain loadings and behavior profiles, which were generated from LOO analysis of CCA, with the weights of the original main mode. We then used the mean and standard deviation of r values to evaluate the reproducibility of CCA in different sample sizes. In CCA1 (ADI&ADOS), the weights of the main CCA mode of each leave-one-out analysis correlated on average above 0.94 with the weights of original main CCA mode in brain loadings and above 0.95 in behavior profiles when the sample was bigger than 122. In CCA2 (SRS&RBS&SSP), the weights of the main CCA mode related on average above 0.92 in brain loadings and above 0.96 in behavior profiles when the sample was bigger than 111. Both CCA analyses are no reproducible for sample sizes smaller than (approximately) 100 subjects. (Figure S2).

Group differences at voxel-wise gray matter volumes
The standard mass-univariate GLM analysis of the VBM data comparing cases and controls did not show significant group differences for voxel-wise GM densities. Figure S3 presents each voxel t-statistics and it is thresholded at uncorrected p<0.05. Figure S3. Results of case-control differences of voxel-wise gray matter densities (p<0.05, uncorrected).

Robustness assessment of the ICA model orders
To assess the reproducibility of components IC10 and IC14 obtained from the mainly reported 100dimensional decomposition, we first correlated the participant loadings from the 100-dimensional factorization with the participant loadings obtained from an alternative factorization. This allowed us to identify the two components from the alternative factorization that are more strongly correlated to IC10 and IC14 respectively, and evaluate then the spatial reproducibility achieved at the alternative factorization. Posthoc statistics, for automatic estimation (91 ICs) and the 50-dimensional IC analysis factorization showed that the composition of components with significant group effects were similar to the original analysis with 100 components. Significant results of both dimensional factorizations are almost equivalent to the IC10 (automatic dimension: p=2.109×10 -4 , 50 dimension: p=0.002) and IC15 (automatic dimension: p=3.557×10 -4 , 50 dimension: p=2.778×10 -4 ) in the analysis of 100-dimensional factorization, however, the ICs corresponding to IC14 in automatic dimension (p=4.733×10 -4 ) and 50 dimensions (p=0.003) did not survive FDR correction (automatic dimension: p<4.733×10 -4 ; 50 dimension factorization: p<0.003) ( Figure S5).   The association analyses were only performed in autism group. ADI, Autism Diagnostic Interview-Revised; RRB, restricted, repetitive behaviors; ADOS, Autism Diagnostic Observational Schedule 2; SRS, Social Responsiveness Scale 2nd Edition; RBS, Repetitive Behavior Scale-Revised; SSP, Short Sensory Profile. a We used SRS parent T-scores. Figure S6. Components with highest loadings in CCA. Panel a shows the three components with highest loadings in CCA1 (correlation with ADI and ADOS subscales). Panel b shows the three components with highest loadings in CCA2 (correlation with SRS, RBS, and SSP). The component maps were thresholded at 3<|Z|<5. Figure S7. The top row shows uncorrected canonical coefficients (uncorrected weights) of the main CCA mode for the CCA1 analyses (ADI&ADOS), and the bottom row for the CCA2 analyses (SRS&RBS&SSP). Panels (a, c) display the degree that each brain component contributed to the main CCA mode in each analysis with respect to the uncorrected canonical coefficients. The two components with significant group effects are displayed in red. Panels (b, d) display the degree that each symptom profile contributes to each analysis Among the uncorrected coefficients, IC14 ranks third among the 100 components when correlating to ADI and ADOS (a), and it ranks fourth in the CCA with SRS, RBS, and SSP (c). CCA, canonical correlation analysis; ADI, Autism Diagnostic Interview-Revised; ADOS, Autism Diagnostic Observational Schedule 2; SRS, Social Responsiveness Scale 2nd Edition; RBS, Repetitive Behavior Scale-Revised; SSP, Short Sensory Profile; IC, independent component.