This study included the acquisition of 662 regionally distributed CT measurements per imaging examination, as extracted by FreeSurfer's recon-all command which processes the input examination with all available atlases. This included extracting measurements of both average and the standard deviation of within-region CT for each supported GM region. This includes all sub-regions of the brain supporting CT measurements across all FreeSurfer supported atlases. Study participants were divided into four groups based on age: early childhood (0-5 years old), late childhood (5-10 years old), early adolescence (10-15 years old), and late adolescence (15-20 years old). We had very few participants greater than 20 years old and so did not include them in a separate group, however, all scatter plots included all participants regardless of age to facilitate visual comparison. We are interested in the diagnostic potential of these clinically acquired measurements and so each measurement (as extracted by FreeSurfer) within each age range was compared in a group-wise manner (autism compared with neurotypical) with receiver operating characteristic (ROC) curve analysis which is summarized with the area under the ROC curve (AUC), Cohen's d statistic (positive/negative values indicate a higher/lower average value in the autistic population relative to the neurotypical population) and a p-value based on the standard t-test for two groups of samples. The p-value was selected as an established method to demonstrate that it is unlikely that our findings were the result of random chance, Cohen's D was selected as it is the most established method to assess effect sizes and the AUC was selected to extend our analysis to the assessment of diagnostic potential. This yielded a total of m = 2,648 group-wise comparisons, yielding a Bonferroni corrected threshold for achieving statistical significance of p < 0.05/m = 1.89e-5.
In order to confirm that the findings reported are the result of group-wise differences between the autistic and neurotypical participants, a statistical model was constructed based on multivariate regression (using the my regress function in MATLAB), adjusting each measurement within each age range in order to control for group-wise differences in age, gender, estimated total intracranial volume and the leading comorbid status of the most common secondary conditions from our two groups: headaches (7% in the autistic group, 19% in the neurotypical group), ADHD (16% in the autistic group, 0% in the neurotypical group), epilepsy (13% in the autistic group, 0% in the neurotypical group), global developmental delay (26% in the autistic group, 0% in the neurotypical group), migraines (3% in the autistic group, 23% in the neurotypical group), and abdominal pain (14% in the autistic group, 11% in the neurotypical group). This model was used to adjust each CT (mean and standard deviation) measurement, in order to evaluate whether group-wise differences between our autistic and typically developing populations are the result of age, gender, intracranial volume or comorbid effects.
Our dataset includes very few exams of participants older than 20 years of age and so this range was not included in the group-wise age-dependent statistical analyses because of the small sample size but was included in all scatter plots for reference. Age-dependent ROC analysis allows us to assess the diagnostic potential of any given FreeSurfer measurement and Cohen's d statistic provides useful effect size measurements at multiple developmental stages. Scatter plots were created to visually present measurements-of-interest from male and female participants as they vary with age. Trend lines in all scatter plots were established with a rolling average (K = 150) implemented in MATLAB (Natick, MA, USA).
Samples of the right insular cortex were as described above. Two age matched frozen postmortem brain tissue specimens, AN13287 and UMB1135, were obtained from National Institute of Child Health Human Development Brain and Tissue Bank for Developmental Disorders (University of Maryland) and Harvard Brain Tissue Resource Center (Boston, MA, USA) through the Autism Tissue Program of Autism Speaks. More detailed information regarding these specimens can be found in our previous work. To control for anatomical variability, landmarks were chosen and marked with India ink during dissections at each brain bank. Tissue blocks (2.0 cm3) were excised from a brain region containing the SVZ region of the lateral ventricles and contained portions of the head of the caudate and corpus callosum from fresh frozen coronal sections. A 5-mm tissue punch containing the SVZ region per sample was used for DNA and RNA extractions. Specimens were de-identified and the study was reviewed by the University of Hawaii (UH) Human Studies Program under application CHS#2016-30171. They determined that this study did not qualify as human subjects research, and thus did not require review and approval by the Human Studies Program or a UH Institutional Review Board (IRB).
DNA was isolated simultaneously with RNA from postmortem SVZ brain samples using AllPrep DNA/RNA/miRNA Universal Kit (80224, QIAGEN) as described above. DNA samples were diluted using EB buffer from the kit to 15 ng/mL and sent to Macrogen, Inc (Seoul, Republic of Korea) for paired-end whole genome Bisulfite sequencing (WGBS). Macrogen performed a quality control check of all samples before preparing libraries by random fragmentation of the DNA and 5' and 3' adapter ligation. Adapter-ligated fragments were PCR amplified and gel purified. Sequencing was performed using the Illumina NovaSeq platform.
Postmortem SVZ human brain tissue samples were thawed on ice and transferred into a 1.5 mL Eppendorf tube. Tissue was homogenized in nuclear extraction buffer (0.32 M Sucrose, 10 mM Tris-HCl pH8.0, 5 mM CaCl2, 3 mM MgCl2, 1 mM DTT, 0.1 mM EDTA and 0.1% Triton X-100) using a disposable pellet pestle (~30 strokes).
Nuclei were isolated by using a discontinuous sucrose gradient. The homogenized tissue solution was transferred into a 50 mL centrifuge tube (3115-0050, Thermo Fisher Scientific) and a sucrose solution (1.8 M sucrose, 10 mM Tris-HCl pH8.0, 3 mM MgCl2 and 1 mM DTT) was added at the bottom to form the gradient. Samples were centrifuged at 75600 x g at 10 C, for 1 h in a Beckman Coulter centrifuge Avanti J-20 XPI with a fixed angle JA-25.50 rotor. Nuclei were stained with a mouse anti-neuronal nuclei (NeuN) Alexa Fluor 488 conjugated monoclonal antibody (Alexa488NeuN antibody) (MAB377X, Millipore) and DAPI (D1306, Thermo Fisher Scientific).
Alexa488NeuN+, DAPI+, and Alexa488NeuNDAPI+ nuclei were sorted using a BD FACSAria III flow cytometer. Parameters and data analysis were established using BD FACSDiva software version 8.0.1 (BD Biosciences, San Jose, CA, USA). Alexa488NeuN antibody fluorescence was detected with the 488 nm laser and a 530/30 filter and DAPI fluorescence was detected using the 407 nm laser and a 450/40 filter. Nuclei were sorted using a 100 mum nozzle at 20 psi with a drop-drive frequency of 29.2 kHz and collected in 1.5 mL Eppendorf tubes with DPBS buffer (14287-080, Gibco).
Postmortem SVZ human brain tissue samples and sorted nuclei were used to extract DNA and RNA using the AllPrep DNA/RNA/miRNA Universal kit (80224, QIAGEN). RNA was treated with DNase and transformed to cDNA using SuperScript  IV Vilo Master Mix Kit (11756050, Thermo Fisher Scientific, Vilnius, Lithuania) according to the protocol of the manufacturer. qPCR assays were performed in a StepOnePlus Real-Time PCR System using 40 ng of cDNA from each sample and TaqMan Fast Advanced Master Mix (4444557, Thermo Fisher Scientific) according to the manufacturer's instructions, using a passive reference of ROX, and analyzed using StepOne Software v2.3 (Thermo Fisher Scientific). TaqMan assays were selected from existing assays available from Thermo Fisher Scientific specific to human with probes spanning exons. Probes labeled as best coverage were preferred. During qPCR, samples were incubated for 2 min at 50 C, then the polymerase was activated at 95 C for 2 min. A total of 40 cycles were performed of 1 s at 95 C and 20 s at 60 C. Once complete, samples were held at 4 C. Co results were normalized to Actinbeta (ActB) amplification, performed in the same assay, and the estimated percent expression was calculated relative to ActB using the delta-Ct method. The TaqMan Assays (Thermo Fisher Scientific) used for this study included: GFAP (Hs00909233_m1), SLC1A3 (Hs00904823_g1), SOX2 (Hs01053049_s1), NES (Hs04187831_g1), ASCL1 (Hs00269932_m1), DCX (Hs00167057_m1), RBFOX3/NeuN (Hs01370653_m1), LIFR (Hs01123581_m1), RBP1 (Hs01011512_g1), FAM107A (Hs00200376_m1), EP300 (Hs00914223_m1), SEPP1 (Hs01032845_m1), DMNT1 (Hs00945875_m1), DMNT3a (Hs00602456_m1), DMNT3b (Hs00171876_m1), and ACTB (Hs01060665_g1).
The WGBS reads were first processed using nf-core/methylseq bioinformatics analysis pipeline to extract the DNA methylation levels for all the samples. In this pipeline, the adaptors and reads with low-quality regions were trimmed or removed by Trim Galore! and the quality of the reads were monitored by FASTQC before and after the trimming. Bisulfite alignment of the reads to hg38 assembly of human reference genome and methylation calling were performed using Bismark. The Bismark coverage reports, containing actual read coverage of detected methylated or unmethylated reads at each position, were used as the inputs for the next step analysis performed by R package RnBeads. The methylation sites with coverage less than 5 reads or with high coverage marked as outliers were removed, while the sites with more than 3 overlapped SNPs or too many missing values were filtered out as well. The rest sites were annotated based genomic regions including CpG sites, tilling regions, genes, and promoters. The differentially methylated region (DMR) analysis can be performed in each level of these region types. Differential methylation on one of the region types was computed based on three metrics: the mean difference in means across all sites in a region of the two groups (unaffected, typically developing "control" vs. ASD) being compared, the mean of quotients in mean methylation, and a combined p-value calculated from all site p-values in the region. A combined rank is computed as the maximum (i.e., worst) value among the three ranks, which were assigned to each region based on the three metrics used. The smaller the combined rank for a region, the more evidence for differential methylation it exhibits. DMRs (genes, promoters) were selected according to the combined rank of a given region and GO enrichment analysis was performed based on these DMRs to predict the functions that might be affected under ASD. Of all sites that were quantified for methylation, we sought to identify sites with significant DNA methylation differences (delta of the beta value) between the ASD and control groups (P < 0.05) using the resampling-based empirical Bayes Methods permutation approach with P < 0.05 and filtered for sites with absolute average methylation differences greater than 10% (delta) between the ASD and control groups. This approach reduces the false discovery rates for non-normally distributed array-based data and offers higher statistical power. Incorporating differences at a threshold of 10% absolute difference in DNA methylation also diminishes the likelihood of random technical errors from true biological differences, increasing the confidence in detecting differences in DNA methylation. We consider these CpGs as differentially methylated loci (DML).
To predict biological functions and gene regulatory networks, differentially methylated or expressed sites identified through WGBS were submitted to the Ingenuity Pathway Analysis software [IPA; QIAGEN; (]. Gene regulatory networks with a log2 fold change of DNAm level were highlighted.