Development of cortical thickness and surface area in autism spectrum disorder.

Autism spectrum disorder (ASD) is a neurodevelopmental disorder often associated with changes in cortical volume. The constituents of cortical volume - cortical thickness and surface area - have separable developmental trajectories and are related to different neurobiological processes. However, little is known about the developmental trajectories of cortical thickness and surface area in ASD. In this magnetic resonance imaging (MRI) study, we used an accelerated longitudinal design to investigate the cortical development in 90 individuals with ASD and 90 typically developing controls, aged 9 to 20 years. We quantified cortical measures using the FreeSurfer software package, and then used linear mixed model analyses to estimate the developmental trajectories for each cortical measure. Our primary finding was that the development of surface area follows a linear trajectory in ASD that differs from typically developing controls. In typical development, we found a decline in cortical surface area between the ages of 9 and 20 that was absent in ASD. We found this pattern in all regions where developmental trajectories for surface area differed between groups. When we applied a more stringent correction that takes the interdependency of measures into account, this effect on cortical surface area retained significance for left banks of superior temporal sulcus, postcentral area, and right supramarginal area. These areas have previously been implicated in ASD and are involved in the interpretation and processing of audiovisual social stimuli and distinction between self and others. Although some differences in cortical volume and thickness were found, none survived the more stringent correction for multiple testing. This study underscores the importance of distinguishing between cortical surface area and thickness in investigating cortical development, and suggests the development of cortical surface area is of importance to ASD.

Introduction

Autism spectrum disorder (ASD) is a spectrum of neurodevelopmental disorders involving deficits in social communication and restricted and repetitive behaviors (American Psychiatric Association, 2013). Magnetic resonance imaging (MRI) studies of brain development in ASD have shown widespread changes in grey and white matter development (Stigler et al., 2011), associated with changes in brain function and connectivity (Amaral et al., 2008, Anagnostou and Taylor, 2011). Recent reviews of the literature have suggested that ASD is associated with subtle increases in brain volume in the first years of life, followed by a period of arrested growth or a decline in volume (Ecker, 2016, Lainhart, 2015, Raznahan et al., 2013). Most of the reviewed studies investigating brain development in ASD focused on cortical volume.

Logically, changes in cortical volume must be due to changes in either cortical thickness or surface area, or both. Cortical thickness is thought to reflect dendritic arborisation and pruning (Huttenlocher, 1990), while surface area may reflect folding and gyrification, both of which depend on division of progenitor cells in the periventricular area during embryogenesis (Chenn and Walsh, 2002, Rakic, 2009). Cortical thickness and surface area have different genetic origins, and distinct developmental trajectories (Panizzon et al., 2009, Raznahan et al., 2011, Wierenga et al., 2014b). As such, it is biologically relevant whether changes in cortical volume are determined by changes in cortical thickness or surface area, as they relate to different underlying biological processes (Ecker et al., 2015).

Studying the development of brain structure is complex as it is associated with regionally specific, non-linear growth patterns and substantial inter-individual variation, even in typical development (Brown and Jernigan, 2012, Lyall et al., 2015, Walhovd et al., 2016, Wierenga et al., 2014b). Furthermore, the literature on brain development varies widely in its methods, including diagnostic criteria utilized, subject characteristics (age, IQ, gender) and MRI acquisition and processing. In addition, ASD is heterogeneous and associated with great variation in symptoms between affected individuals. This complicates the comparison of findings between studies and is a likely cause for contradictory findings so far (Gronenschild et al., 2012, Raznahan et al., 2013, Stigler et al., 2011, Walhovd et al., 2016). Longitudinal studies of a reasonable size with well-matched control groups can help address these issues.

Some studies have already investigated developmental changes in cortical thickness or surface area longitudinally. Hazlett and colleagues reported greater cortical volume in ASD, with similar thickness, and suggested that surface area may therefore be driving the volume increases in young children with ASD, although they did not assess it directly (Hazlett et al., 2011). However, a recent study assessed both cortical thickness and surface area in similarly aged children (2–9 years old) and found no differences in surface area, but rather widespread differences in thickness (Smith et al., 2016). Here, thickness decreased in typically developing children, but not in children with ASD. A study of older children and adolescents (8–15 years old) found steeper developmental decreases in thickness in ASD, primarily in temporal and occipital regions (Hardan et al., 2009). Another study reported that greater developmental decreases in thickness continued in ASD between the ages of 14 and 24 in left ventral temporal cortex and superior parietal cortex, while the development of surface area did not differ from controls (Wallace et al., 2015).

Zielinski and colleagues investigated cortical thickness longitudinally, in a large sample of subjects with ASD with a broad age range (3–39 years old). They reported increased thickness in early childhood, followed by accelerated thinning into later childhood and adolescence and decelerated thinning in early adulthood. In adolescence, reduced thickness was found predominantly in frontal areas, with a more protracted timeline for more posterior regions (Zielinski et al., 2014).

One of the methodological challenges in the field of developmental imaging is that brain regions are not independent, which makes estimating the true degrees of freedom for conventional multiple comparison corrections inaccurate. There is an ongoing debate on how best to address this issue. Highly stringent corrections such as Bonferroni are overly conservative, as they assume the measures under investigation are independent. In this study, we applied a fairly lenient correction, typical for the neuroimaging literature (with an increased risk of type II errors). In addition, we borrowed a singular value decomposition-based correction method for multiple comparisons from the genetics field that takes mutual dependence of the measures of interest into account (Nyholt, 2004).

In conclusion, the studies above suggest that cortical development differs in ASD from controls throughout childhood and adolescence, where differences are mainly ascribed to cortical thickness. However, the number of studies investigating cortical thickness and surface area simultaneously is still small (Smith et al., 2016, Wallace et al., 2015), as is the number of subjects these studies were able to include, and none have yet used a longitudinal design. Therefore, we set out to investigate the developmental trajectories of cortical thickness and surface area in a sample of 180 subjects from late childhood to young adulthood (9–20 years old), using an accelerated longitudinal design.

Methods and materials

The Institutional Review Board of the University Medical Center Utrecht, the Netherlands, approved the study protocol. Subjects aged 18 years or older provided written informed consent after full disclosure of the study purpose and procedure. For subjects under the age of 18 years, the parents provided written informed consent and subjects gave assent

Participants and clinical measures

We acquired 406 whole-brain MRI scans from 286 subjects (132 meeting DSM-IV (American Psychiatric Association, 2000) criteria for autism or Asperger syndrome and 154 typically developing control subjects). To account for group differences in age and gender we matched individuals on these factors using parametric models which are provided by the R package “MatchIt”. This matching procedure resulted in a sample of 180 subjects (90 in each group), aged between 9 and 20 years old (Table 1). The ASD group included a total of 115 scans, including 24 individuals with two or more scans. The control group had a total of 141 scans, with 51 individuals were scanned twice or more.

Table 1

Demographic and clinical characteristics of the participants.

		ASD	Controls	Group
		(N = 90)	(N = 90)	Differences
Gender	Female/male	16/74	16/74	n.s.
Age at scan	N scans wave 1	90	90
	Years M (SD)	13.8 (2.9)	13.2 (3.0)	n.s.
	N scans wave 2	22	36
	Years M (SD)	15.3 (2.4)	15.0 (2.9)	n.s.
	N scans wave 3 +	3	15
	Years M (SD)	15.0 (1.2)	16.1 (1.7)	n.s.
Hand preference	N right-handed/other	72/18	70/20	n.s.
SES	Paternal education years M (SD)	13.6 (2.6)	14.0 (2.2)	n.s.
	Maternal education years M(SD)	13.1 (2.4)	13.8 (2.2)	p = 0.04
Total IQ	M (SD) at baseline	103.7 (19.1)	114.1 (15.2)	p = 0.0001
ADI	Social deficits M (SD)	18.8 (4.8)	–
	Communication deficits M (SD)	14.6 (4.2)	–
	Ritualistic-repetitive behavior M (SD)	4.9, (2.5)	–

ASD; autism spectrum disorder; N, number; n.s., not significant; M, mean; SD, standard deviation; SES, socio-economic status; IQ, intelligence quotient; ADI, autism diagnostic interview

Diagnosis for autism or Asperger syndrome was clinically established by a psychiatrist, and was confirmed by trained and qualified clinicians using the Autism Diagnostic Interview-Revised (ADI-R) (Le Couteur et al., 2003). To confirm absence of psychiatric illness in the control group, we administered the parent version of the Diagnostic Interview Schedule for Children (DISC, version 2.3 or IV) (Shaffer et al., 2000) for individuals under 18 years of age. Individuals 18 years or older were administered the Mini-International Neuropsychiatric Interview (MINI) (Lecrubier et al., 1997).

In addition, total IQ was estimated using a short Dutch version of the Wechsler intelligence scales (WISC-R/WISC-III or WAIS-III dependent on the age of the participant), including the subtests Vocabulary, Block Design, Similarities and Object Assembly.

Controls were excluded in case of psychiatric morbidity or first-degree relatives with a history of psychiatric problems. In both groups, additional exclusion criteria were IQ below 70, any major physical or neurological illnesses or the presence of metals in the body that precluded the MRI session. Table 1 lists participant characteristics for both samples.

Prior to the MRI scan, children under 13 years of age were acclimated to the MRI procedure in a practice session using a dummy scanner as described previously (Durston et al., 2009); adolescents over 13 years were also offered the opportunity to participate in such a practice session.

MRI acquisition

Scans were acquired using two identical 1.5-T MRI-scanners (Philips, Best, The Netherlands). Scanning sessions were randomly assigned to one of the scanners, post hoc tests showed no effects of scanners for any of the ROIs. A T1-weighted three-dimensional fast field echo scan of the whole head was acquired with 130 to 150 1.2-mm contiguous coronal slices (echo time 4.6 ms; repetition time 30 ms; flip angle 30°; field of view 256 mm; in-plane voxel size 1 mm × 1 mm). Independent neuroradiologists evaluated all MRI scans and no gross morphological abnormalities were reported for any of the participants.

MRI processing

The T1-weighted images were processed using the well-validated open access software package FreeSurfer v5.1.0 (Fischl, 2012) (available at http://surfer.nmr.mgh.harvard.edu/). This software package involves an automated procedure that allows quantitative assessment of brain anatomy including volumetry of subcortical structures and cortical morphometry, with accuracy comparable to manual methods (Fischl et al., 2002). The procedures have been described previously (Dale et al., 1999, Fischl and Dale, 2000, Fischl et al., 2004a, Fischl et al., 2004b, Fischl et al., 2002, Fischl et al., 1999). In short, brain segmentation consists of removing non-brain tissue using a deformable template model and neuroanatomical labeling based on both voxel intensity values and a probabilistic atlas (Fischl et al., 2002). Next the cortical surface is reconstructed where the segmentation of white matter is used to derive a surface representing the grey-white matter boundary (Dale et al., 1999). Finally, by incorporating both geometric information derived from the cortical model and standard neuroanatomical conventions (Desikan et al., 2006, Destrieux et al., 2010), the procedure automatically assigns a neuroanatomical label to each location on the cortical surface. All scans were processed on the same hardware using the same software to avoid potential confounds (Gronenschild et al., 2012).

All scans were processed using the longitudinal processing pipeline of FreeSurfer (Reuter and Fischl, 2011). This method reduces the risk of potential over-regularization of longitudinal image processing (Reuter et al., 2012).

The output of each scan was checked by two independent raters (VM/LW) who were blind to subject identity or group assignment. Errors were corrected where necessary following the standardized procedures documented on the FreeSurfer website (http://surfer.nmr.mgh.harvard.edu/). The types of errors that frequently required user editing were incomplete skull stripping and mis-classification of the white matter segmentation.

For the study of cortical morphometry, we analyzed 34 areas per hemisphere mapped with the Desikan-Killiany atlas (Desikan et al., 2006). For each of these regions we assessed cortical volume (mm3), surface area (mm2) and thickness (mm). For each area, surface area was measured along the white surface; volume and thickness were measured as the volume and the average distance, respectively, between the parcellated portions of white and pial surfaces (Fischl and Dale, 2000). Total, left (lh) and right (rh) hemisphere values were obtained for cortical volume, surface area (CS) and thickness (CT) by summing or averaging each measure across all areas included. Average thickness throughout the cortex was computed applying the formula: Total CT = [(lh.CT ∗ lh.CS) + (rh.CT ∗ rh.CS)] / lh.CS + rh.CS.

Statistical analysis

The developmental trajectories of each measure of interest were modeled with the Linear Mixed Model as implemented in the lme4 library (Bates et al., 2012) in the R statistical package (R Core Team, 2013). This method allows for inclusion of multiple time points per subject, while accounting for unbalanced data structure of irregular time intervals between scans and unequal numbers of scans per subject, e.g. differences in number of follow up scans between children with ASD and controls.

The best-fit model was selected in two phases based on previous methods as previously described (Wierenga et al., 2014a, Wierenga et al., 2014b). First, the best fit growth model (cubic, quadratic, linear) was selected for all participants, including both children with ASD and controls. We used a step-down selection procedure, starting with the full model including all age terms. Each brain measure of interest of the ith individual at the jth time point was modeled as follows:where dij represents the within-person dependence and the eij term is residual error. Age effects were fixed, while the intercept and the dij term were modeled as random effects. If the cubic age effect was not significant at p < 0.05, it was removed from the model and we next tested the model including linear and quadratic age terms only, and so on. In the next step, we examined whether the growth model differed between children with ASD and controls, by including main and group interaction terms. For example, if the cubic model was selected in step 1 the group model was modeled as follows:

Again the selection procedure started with the full model including both group and interaction terms. The best-fit model was selected when the estimated interaction coefficient of the highest order (in this example β7) would reach a threshold of p < 0.05. If not, we stepped down to the model including main effects only, and tested whether the estimated main effect of group was significant (p < 0.05). If this was not the case we concluded there were no significant differences between the groups. By using a Bayesian simulation we estimated that our sample size has a probative power of at least 80% to detect a 5% difference in brain volume between children with autism and controls at a false discovery rate of 5%.

In addition, we calculated a more stringent correction for multiple comparisons. We used a singular value decomposition-based correction method that is often applied to correct for multiple testing of single-nucleotide polymorphisms (Nyholt, 2004) and that has shown to be applicable in neuroimaging studies before (Hedman et al., 2015, Mandl et al., 2013). This method applies a correction factor on the degrees of freedom by taking the dependencies between brain regions into account. Nyholt et al. showed that by measuring the variance of eigenvalues of the correlation matrix of raw measurements the effective number of degrees of freedom can be estimated. For example, if the measurements are completely independent then all eigenvalues are equal to one and the variance is zero. In the case that the measurements are completely dependent then the first eigenvalue equals the number of variables in the correlation matrix and the rest of the eigenvalues is zero (and hence the variance is maximal). The correction factor was computed separately for cortical volume, thickness, and surface area values. Next, the average correction factor (59) was used to calculate the corrected degrees of freedom.

Results

Total cortical volume, surface area and average thickness

The best fitting growth models for total cortical volume, surface area and thickness are shown in Fig. 1 and the regression coefficients are reported in Table 2. For total cortical volume, the best fit was a cubic developmental trajectory, with stable, nearly 3% smaller volume in subjects with ASD than typically developing controls. For surface area, we found an age ∗ group interaction: children with ASD had almost 5% smaller total surface area at younger ages, but comparable surface area by late adolescence. Average cortical thickness linearly declined with age for both groups, with no between-group differences.

Fig. 1

Developmental trajectories for total cortical volume (in mm3), total cortical surface area (in mm2), and average overall cortical thickness (in mm) for subjects with ASD (green) and controls (grey). Age (in years) is depicted on the x-axis.

Table 2

Global and regional differences between groups.

Total values	Model	Effect	Coefficient (s.e.)	p
Total cortical volume	Cubic	ASD < controls	− 16,168 (7351.007)	0.029⁎
Total cortical surface area	Linear	Interaction	− 4058 (2274.939)	0.006⁎⁎
Average cortical thickness	Linear			n.s.
Total subcortical grey volume	Quadratic	ASD < controls	(0.000)	0.013⁎

Regression coefficients and p-values of main group and group ∗ age interaction terms are provided in the table. ASD, autism spectrum disorder; TD, typical development; LH, left hemisphere; RH, right hemisphere; n.s., not significant.

p < 0.05.

p < 0.01.

p < corrected for multiple comparisons using singular value decomposition

Regional differences in cortical volume, surface and thickness

We found differences in cortical volume between children with ASD and their typically developing peers in 21 of 68 cortical regions of interest (p < 0.05, uncorrected). We found main effects of group on cortical volume in left fusiform gyrus, bilateral inferior parietal gyrus, left inferior temporal gyrus, bilateral middle temporal gyrus, right precentral gyrus, left rostral anterior cingulate, bilateral superior frontal gyrus, bilateral marginal gyrus, bilateral insula and right precuneus. In these regions, cortical volume was decreased by 3 to 5% in children with ASD. The regression coefficients for main group effects and group ∗ age interaction effects are reported in Table 2. To test the robustness of these findings, we also fit the model for each group separately and ran gamm spline fitting analyses. The results were highly similar to those reported here.

For 15 of these 21 regions, between-group differences in volume were associated with differences in cortical surface area rather than cortical thickness, with the exception of right and left inferior temporal gyrus and right posterior cingulate gyrus. These latter regions showed greater differences between the developmental trajectories of cortical thickness than surface area between children with ASD and controls. There were four regions where we did not find group differences in cortical thickness or surface area, despite between-group differences in volume. Regional differences in surface area were related to either main group effects, where children with ASD had smaller surface area than controls, or group ∗ age interaction effects, where children with ASD showed a smaller decline with development than their typically developing peers or even developmental increases in surface area. These regional differences are shown in Fig. 2.

Fig. 2

Regional effects of group (blue) and trajectory differences between groups (green and purple) across the cortical mantle. Effects are displayed for the left and right hemispheres in cortical volume (top row) cortical surface area (middle row) and average cortical thickness (bottom row). Regions that showed smaller values in ASD compared to controls across the entire age range are indicated in blue. Area's that showed group ∗ age interaction effects are displayed in green and purple. Most of these regions showed linear age effects with smaller declines in the ASD group compared to typically developing peers. Exceptions are the volumes of the left superior parietal gyrus and right lingual gyrus as well as the thickness of the left superior frontal gyrus that all showed cubic, rather than linear age effects.

In addition to these exploratory analyses at a lenient statistical threshold, we also investigated which of these results survived correction for multiple comparisons using singular value decomposition; only three results remained, all reflecting between-group differences in the developmental trajectories of surface area: in the left banks of superior temporal sulcus, left post central gyrus and right supramarginal gyrus. These regions all showed group ∗ age interaction effects, with a smaller decline with development in ASD than in typical development (see Fig. 3).

Fig. 3

Developmental differences in surface area (in mm2) for subjects with ASD (green) and controls (grey) in three regions that survive stringent correction for multiple comparisons. Age (in years) is depicted on the x-axis.

Discussion

We report that differences in cortical development in ASD are largely driven by differences in the development of cortical surface area. We found that in typical development, surface area decreased from late childhood to early adulthood, whereas it stayed relatively stable in ASD, both globally, and in most cortical regions where we found between-group differences. The only between-group differences that survived stringent correction for multiple comparisons were seen in surface area, no differences remained significant in volume or cortical thickness.

Our findings add to other, large longitudinal studies that have investigated the development of cortical volume (Lange et al., 2015) and thickness (Zielinski et al., 2014) in ASD. These studies, and two more recent, smaller studies investigating cortical surface area and thickness (Smith et al., 2016, Wallace et al., 2015), have suggested changes in the development of cortical thickness in ASD, but not surface area. In contrast, we report that changes in the development of cortical volume are more likely related to changes in the development of surface area. The age-range of subjects included in these studies varied widely, with Smith and colleagues including only subjects that were much younger than those in the current paper (Smith et al., 2016). Wallace and colleagues included a sample that overlapped with the current one in terms of the age-range included. Yet, the youngest age in their sample was 14 years, in mid-adolescence, which of course may have impacted which statistical model gave the best fit.

Most of the regional between-group differences in volume we found were associated with differences in the development of surface area. When we used a stringent correction for multiple testing that takes the interdependence of regions into account, significance was maintained for three regions. All of them showed the pattern of developmental differences in surface area that was most commonly found in other regions at less stringent statistical thresholds. The three regions that survived stringent multiple comparison correction were the left banks of superior temporal sulcus, postcentral area, and right supramarginal area. These regions have all been implicated in social functioning: in a meta-analysis of functional MRI studies in ASD, left postcentral region was associated with complex social cognition in children and adolescents with ASD and left banks of superior temporal sulcus area were associated with auditory perception and interpretation of social stimuli (Philip et al., 2012). In a large, recent study on resting state functional connectivity, these same areas were also associated with face and motion processing related to social behavior (banks of the superior temporal sulcus) and somatosensory functioning and sense of self (postcentral gyrus) (Cheng et al., 2015). Right supramarginal area has been associated with the distinction between self and others (Hoffmann et al., 2015).

When we applied a more liberal statistical correction, more similar to other neuroimaging studies, we found more differences between groups. These included differences in cortical volume and cortical thickness, which is more in line with the studies mentioned above.

We found between-group differences in volume for 21 regions. For the majority of these (16/21), the developmental trajectories for cortical volume were parallel for typically developing controls and subjects with ASD, with a lower volume associated with ASD. This is in keeping with earlier findings of reduced cortical volume in ASD during adolescence (Lange et al., 2015). In the other five regions, we found either linear age ∗ group interactions, with controls showing a steeper decline in volume, or a cubic pattern. These findings were located in regions associated with changes in social cognition, audiovisual processing, and executive and motor functioning in ASD (Philip et al., 2012). Furthermore, we found differences in the development of surface area in 15 of these 21 regions. In five of these regions, ASD was associated with stable reductions in surface area compared to controls. However, for the majority of these regions (10/15), we found an age ∗ group interaction, with typically developing controls showing steeper decline than subjects with ASD. This is the pattern that we also observed in the three regions that survived the more stringent correction for multiple comparisons. This suggests that the pattern of change we found in the three regions that survived the more stringent statistical correction may be representative of a more widespread pattern. We found differences in the development of cortical thickness in three of 21 regions, but only when using the more liberal correction for multiple comparisons. This suggests that these findings are more likely to be erroneous. These results highlight the importance of statistically controlling for multiple comparisons in this and similar studies.

In interpreting our results, some limitations should be considered. Our scans were acquired at 1.5 T, where many recent studies have used 3 T scanners. In spite of similar studies also employing a 1.5 T field strength (Smith et al., 2016) and studies showing results from FreeSurfer are comparable at 1.5 T and 3 T, a higher field strength would inevitably have increased signal to noise ratio. This is particularly relevant to cortical thickness, as this measure has been suggested to be particularly sensitive to SNR (Han et al., 2006). Furthermore, even larger samples than the one included here, with more longitudinal scans per participant, and covering a broader age range, are necessary to more accurately estimate the developmental trajectories of cortical volume, thickness and surface area. Finally, our ASD subjects were high-functioning, limiting generalizability to lower-functioning individuals in the broader autism spectrum.

In conclusion, we found that the main difference in cortical development between typical development and ASD was in surface area, where a pattern of decline in typically developing controls between the ages of 9 and 20 was absent in ASD. This pattern survived stringent correction for multiple comparisons in the left banks of superior temporal sulcus, postcentral area, and right supramarginal area.

Conflict of interest

All authors declare that they do not have any conflicts of interest.