Broad transcriptomic dysregulation across the cerebral cortex in ASD

Classically, psychiatric disorders have been considered to lack defining pathology, but recent work has demonstrated consistent disruption at the molecular level, characterized by transcriptomic and epigenetic alterations.1–3 In ASD, upregulation of microglial, astrocyte, and immune signaling genes, downregulation of specific synaptic genes, and attenuation of regional gene expression differences are observed.1,2,4–6 However, whether these changes are limited to the cortical association areas profiled is unknown. Here, we perform RNA-sequencing (RNA-seq) on 725 brain samples spanning 11 distinct cortical areas in 112 ASD cases and neurotypical controls. We identify substantially more genes and isoforms that differentiate ASD from controls than previously observed. These alterations are pervasive and cortex-wide, but vary in magnitude across regions, roughly showing an anterior to posterior gradient, with the strongest signal in visual cortex, followed by parietal cortex and the temporal lobe. We find a notable enrichment of ASD genetic risk variants among cortex-wide downregulated synaptic plasticity genes and upregulated protein folding gene isoforms. Finally, using snRNA-seq, we determine that regional variation in the magnitude of transcriptomic dysregulation reflects changes in cellular proportion and cell-type-specific gene expression, particularly impacting L3/4 excitatory neurons. These results highlight widespread, genetically-driven neuronal dysfunction as a major component of ASD pathology in the cerebral cortex, extending beyond association cortices to involve primary sensory regions.

Broad attenuation of transcriptomic regional identity 2 3 We previously observed an attenuation of typical gene expression differences between two regions, 4 frontal and temporal lobe in ASD, 4,5 which we refer to here as an "Attenuation of Transcriptomic Regional Identity" 5 (ARI). To assess whether this was a broader phenomenon, we systematically contrasted all unique pairs of 11 6 cortical regions (55 comparisons in all) using a conservative statistical approach to account for differences in 7 sample size across regions, while correcting stringently for multiple comparisons (Fig. 2a, Methods). We further 8 validated the identified transcriptomic regional identity patterns in our control samples with those from an external 9 data source, the Allen Brain Atlas 9 (Supplementary Table 4, Extended Data Fig. 5, Methods). Ten pairs of 10 regions exhibited significantly greater ARI patterns in ASD compared to controls, with an additional 31 out of the 11 55 pairs of regions exhibiting a trend towards attenuation in ASD (Fig. 2b, Supplementary Table 4, Extended 12 Data Fig. 5, Methods). These results provide evidence in support of widespread ARI across the cerebral cortex 13 in ASD for the first time, across both gene and isoform levels (Extended Data Fig. 5). Additionally, we observed 14 a regional anterior -posterior gradient, with nine of the ten region pairs exhibiting significant ARI in ASD 15 containing either BA17 or BA39-40 ( Fig. 2c-d).
Notably, BA17 was also one of the regions with the largest case-16 control differences in gene expression. To determine how gene expression changes were dispersed across 17 regions in these pairs, we used a conservative filtering process to identify individual genes exhibiting ARI 18 (Methods, Supplementary Table 4). Although these genes were widely dysregulated, the posterior regions 19 BA17 and BA39-40 exhibited the greatest changes ( Fig. 2c-d, Extended Data Fig. 6). ARI genes were also 20 comparably disrupted in the dup15q samples (Extended Data Fig. 6), suggesting that transcriptomic regional 21 identity attenuation in the cerebral cortex is shared across heterogenous forms of ASD. 22 To identify the biological processes contributing to ARI gene dysregulation in ASD, we grouped together 23 all of the ARI genes that were either downregulated (1,881 genes) or upregulated (1,695 genes) with a 24 pronounced posterior effect in ASD (Methods). The downregulated set of ARI genes showed broad enrichment 25 for neuronal cell-type-specific markers and RNA processing pathways, and contained many transcription factors 26 (Fig. 2c, Supplementary Table 4). The upregulated ARI genes also contained many transcription factors and 27 were enriched for oligodendrocyte progenitor cell (OPC) and astrocyte cell-type markers along with metabolic 28 and development pathways. ARI gene dysregulation was further characterized by subsequent co-expression 29 network analysis, which further refined the topology and pathways involved. 30 31 Refining disrupted gene co-expression networks in ASD 32 33 We next used weighted gene correlation network analysis (WGCNA) 10 across all samples to partition 34 genes into co-expression modules capturing potentially shared biological functions or regulation (Methods). We 35 identified a total of 35 gene modules, of which 9 were downregulated and 15 were upregulated in ASD 36 (Supplementary Table 5-6, Extended Data Fig. 7). We further generated networks using isoform-level 37 quantifications, identifying 61 isoform modules. Of these, 39 were distinct from the gene modules, with 5 38 downregulated and 9 upregulated in ASD (Supplementary Table 5-6, Extended Data Fig. 8). In total, 38 gene 39 and isoform modules were dysregulated in at least one region in ASD. Most of these fell into two broad groups 40 -either dysregulated (1) cortex-wide with comparable magnitude across regions, or (2) with significantly variable 41 magnitude across regions. Again, dup15q effects were similar to ASD effects, but were greater in magnitude 42 (Extended Data Fig. 7-8, Supplementary Table 6).  Table 6). These include GeneM9, an upregulated neuronal module with a significant enrichment for non-coding 47 5 genes; GeneM32, a strongly upregulated reactive astrocyte module with the greatest overall magnitude of 1 dysregulation; and GeneM24, a downregulated module enriched for endothelial and pericyte marker genes which 2 are involved in blood-brain-barrier functions (Fig. 3b, Extended Data Fig. 7, Supplementary Table 6). These 3 modules replicate previous findings of neuronal upregulation, astrocyte reactivity, and BBB disruption in ASD, 1,4-4 6 but extend these findings by demonstrating that these processes are widespread across the cerebral cortex. 5 Two modules -GeneM5 and IsoM37 -demonstrated cortex-wide dysregulation along with significant 6 enrichment for ASD-associated common genetic variation (Fig. 3b-d). 11 GeneM5 is down-regulated in ASD, 7 contains many neuronal genes involved in synaptic plasticity, and significantly overlaps with the synaptic module 8 CTX.M16 previously identified by Parikshak et al. 5 (Fig. 3a, Fig. 3d, Supplementary Table 5-6). In addiiton to 9 common genetic variation, GeneM5 is also significantly enriched for genes containing rare de novo protein 10 disrupting mutations associated with ASD, including the high-confidence risk genes GRIN2A, MYO5A, and 11 BTRC 12 (Supplementary Table 5-6, Methods). This demonstrates another point of convergence of rare and 12 common risk variants on shared biological processes in ASD. 39 GeneM5 is enriched in cortical lower layer 4-6 13 excitatory neuron cell-type markers (Extended Data Fig. 7), 13 identifying them as a point of convergence for 14 rare and common genetic risk in ASD. Finally, IsoM37 is enriched for ASD common genetic risk variants (but not 15 rare mutations), is upregulated in ASD, and contains genes involved in protein folding (Fig. 3a, Fig. 3c, 16 Supplementary Table 6). To our knowledge, this is the first report of an upregulated ASD transcriptomic 17 signature that is associated with known ASD risk variants. 18 Magnitude of effect parallels anterior-posterior gradients 19 In addition to observing profound cortex-wide dysregulation in ASD, we found 13 modules that exhibited 20 their most pronounced ASD effect in BA17, as measured against a permuted distribution containing all regions 21 (Extended Data Fig. 7, Supplementary Table 6, Methods). Of these, 12 showed significant enrichment for ARI 22 genes (half up-regulated and half down-regulated in ASD) and all 13 had anterior -posterior gradients of 23 expression in neurotypical samples, indicating that these modules contribute to transcriptomic regional identities 24 that are observed in neurotypical controls, but attenuated in ASD (Extended Data Fig. 7, Supplementary Table  25 6). Six of these modules were more highly expressed in posterior regions in neurotypical subjects and were 26 observed to be downregulated in ASD across the cortex (Fig. 4a, Extended Data Fig. 7, Supplementary Table  27 6). These include GeneM23, an oligodendrocyte-specific module consisting of genes important for organelle 28 regulation and intracellular restructuring; GeneM14, a neuronal module that contains genes involved in neurite 29 morphogenesis and is also strongly downregulated in BA41-42-22; and GeneM3, a neuronal module enriched  30 for energy generation and neuronal processes that are highly energy dependent, such as vesicle transport ( Fig.  31 4b-c, Supplementary Table 5-6). GeneM3 is also significantly enriched for cell-type markers specific to layer 32 4-5 excitatory neurons (Extended Data Fig. 7). 13 The next four modules were more highly expressed in anterior 33 regions in neurotypical subjects and exhibited cortex-wide upregulation in ASD (Fig. 4a, Extended Data Fig. 7, 34 Supplementary Table 6). These include GeneM8, a microglial module containing genes involved in immune 35 signaling and phagocytosis; and GeneM7, an immune response module containing genes such as NF-kB and 36 interferon response pathways (Fig. 4b, Supplementary Table 5-6). Although neuronal and oligodendrocyte 37 downregulation along with immune and microglia upregulation have been previously reported in ASD, 1,4-6 these 38 findings indicate that this dysregulation is widespread across the cerebral cortex, with increased magnitude in 39 posterior regions, a pattern most pronounced in BA17. 40 The last three modules were observed to be significantly upregulated in ASD only in BA17. One of these, 41 GeneM4, is an inhibitory neuron module containing many genes important for various intracellular signaling and 42 maturation processes, such as SCN9A (Fig. 4a- Fig. 7). We also identified four other modules exhibiting strong region-specific 46 dysregulation in regions other than BA17 (Extended Data Fig. 7). For example, the module GeneM34, which 47 contains genes involved in cellular stress response regulatory processes, is upregulated with the greatest 1 magnitude in BA4-6 and shows no significant effect in BA17 (Extended Data Fig. 7, Supplementary  None of the gene modules with regionally-variable magnitudes of ASD effect were significantly enriched for  3  known ASD genetic risk variants.  4 Cell-type changes mirror regional variation 5 We finally sought to determine what might be driving the observed changes in magnitude of ASD effect 6 across regions. It is well established that BA17 is the most neuronally dense region in the human brain, with a 7 notable expansion in the thickness of L3/4, compared with other cortical regions. 18 Likewise, there is an anterior-8 posterior gradient in neuronal density observed in mice and primates. 14-17 As such, we posited that regional 9 variation in cell density could be contributing to regional differences in magnitude of ASD effect. Regional 10 neuronal density across multiple brain regions has not been quantitatively studied in the human brain, but such 11 gradients have been established across some regions in non-human primates. 15,16 Therefore, we compared the 12 region-specific ASD effect size changes in our gene modules to regional neuronal nuclei density measured in 13 primates 15 for 6 matched regions across species. We observed a significant association between neuronal 14 density and the effect sizes for several modules dysregulated in ASD (seven with FDR < 0.05, and an additional 15 eight with FDR < 0.1, Extended Data Fig. 9, Supplementary Table 7). Further, L3/4 thickness was also 16 associated with the region-specific ASD effect sizes in dysregulated modules (Supplementary Table 7). 17 These observations motivated us to perform single-nucleus RNA sequencing (snRNA-seq) in a small 18 cohort of individuals to help evaluate how distinct neural cell-types could be contributing to the regional variance 19 in ASD transcriptomic dysregulation identified with bulk RNA-seq (Fig. 4d, Extended Data Fig. 9, 20 Supplementary Table 7, Methods). We sequenced over 150,000 nuclei from ASD and control samples across 21 frontal and occipital cortices with matching bulk RNA-seq. From these data, we identified 35 distinct cell clusters 22 and 4,953 cell-type-specific DE genes in ASD subjects in the frontal and occipital cortex. The vast majority of 23 these were DE in excitatory neurons in both regions, and exhibited larger effects overall in the occipital lobe ( Fig.  24 4f). While statistical power limited our ability to detect significant cell-type proportion differences between regions 25 or diagnoses (Methods), we do observe that excitatory neurons are increased in proportion by ~5% in BA17 26 across both control and ASD subjects compared to frontal regions (Extended Data Fig. 9), corresponding with 27 the primate neuronal density measurements. To predict how cell-type proportions may vary across our entire 28 bulk RNA-seq dataset, next we utilized cell-type markers from our snRNA-seq to perform cell-type deconvolution 29 in all samples (Methods). We identified 11 significant cell subtype proportion changes present across six 30 different regions in ASD, characterized by neuronal decreases and astrocyte and microglia increases (Fig. 4e,  31 Extended Data Fig. 9, Supplementary Table 7). We also found many anterior-posterior cell-type proportion 32 gradients in control subjects that are attenuated in ASD (Extended Data Fig. 9, Supplementary Table 7), 33 mirroring patterns observed with our bulk RNA-seq transcriptomic regional identity analysis. 34 When directly comparing cell-type-specific DE and deconvolved proportional changes in ASD with 35 regional variability in the larger bulk transcriptome sample, we observed a convergent signal within excitatory 36 neuronsin particular, those in L3/4 (ExNeuron4; Fig. 4e-f, Extended Data Fig. 9). Recapitulating the known 37 increase in thickness of L3/4 in BA17 compared with other cortical regions 18 , we observed a significant increase 38 in the estimated proportion of the ExNeuron4_L3/4 cluster in posterior regions, peaking in BA17 ( Fig. 4e;  39 Extended Data Fig. 9g). This regional pattern was significantly attenuated in ASD, with a ~2-fold median 40 reduction in estimated Ex4 neuronal proportion in BA17 compared with controls. This cell cluster, with marker 41 genes RORB, PCP4, CUX2, PHACTR2, and EYA4, also exhibited substantially greater cell-type-specific DE in 42 snRNA-seq profiling from BA17 compared with frontal cortex (90 vs 0 DE genes, respectively; Extended Data 43 Fig. 9d). Similarly, BA17 shows a substantially greater upregulation of inhibitory neuron genes in the single cell 44 data, consistent with the observed greater up-regulation of GeneM4 (inhibitory neuron) in BA17 (Fig. 4f). 45 Substantial changes in gene expression are also evident in other cell subtypes (Extended Data Fig. 9, 46 Supplementary Table 7), such as microglia_2, which shows a strong and specific increase in DE genes in ASD 47 7 BA17 compared to frontal regions. These observed intracellular/cell-type changes in neuronal and microglial 1 gene expression are further supported by another snRNA-seq dataset containing a small ASD cohort, which 2 assessed a single region. 19 Here, through performing multi-region snRNA-seq and cell-type deconvolution, we 3 show that predicted cell-type proportions as well as cell-type-specific gene expression profiles are impacted 4 across the ASD cerebral cortex. Importantly, we see increased cell-type-specific transcriptomic dysregulation 5 and lowered neuronal proportions with a notable convergence within L3/4 excitatory neurons in ASD BA17, a 6 region where neuronal proportions are neurotypically abundant. These changes likely contribute to the 7 pronounced ASD effect we observe with bulk RNA-seq in this region. 8 9 Discussion 10 11 Overall, the findings presented here substantially expand our understanding of ASD pathology beyond 12 the previously established 'downregulated neuron' and 'upregulated glia/immune' functional categories observed 13 in frontal and temporal lobe. We identify gene and isoform expression changes in ASD that extend across the 14 cerebral cortex, many neural cell-types, and specific biological processes (Extended Data Fig. 10), including 15 primary sensory areas in addition to association areas. 1,[4][5][6] We find that the recently observed reactive astrocyte 16 upregulation and blood-brain barrier membrane transport downregulation 1 is extended cortex-wide in ASD. It is 17 interesting to speculate that the substantial changes observed in area 17, a primary sensory region, may be 18 related to the widespread observation of sensory hypersensitivity or processing abnormalities in ASD. 40 19 Nevertheless, this region shows the most profound changes in gene expression, a clear demonstration that these 20 alterations are not specific to higher association areas. 21 Furthermore, we find that other dysregulated pathways observed before in ASD -particularly upregulated 22 immune response and reactive microglia genes, along with downregulated neurite morphogenesis and neuronal 23 energy pathway genes -are not only impacted cortex-wide in ASD, but impacted in a regional gradient that 24 reflects fundamental elements of cortical cytoarchitecture, such as neuronal density. It is also notable that the 25 magnitude of region-level differences in ASD parallels regional variance in attenuation of transcriptomic identity 26 (ARI gene dysregulation), suggesting that they reflect related processes. That the gradient of region-specific 27 changes between ASD and controls coincides with both neuronal proportion differences and cell-type-specific 28 transcriptomic dysregulation further suggests that the interplay of cytoarchitecture and cell-type gene expression, 29 rather than a single one of these features, influences our ability to observe transcriptomic changes in bulk tissue. 30 Given the connection between regional cytoarchitecture, local circuits and long-range brain connectivity, 20,21 31 parsimony suggests that in addition to developmental patterning contributions, 5,22 the diminution of transcriptomic 32 regional identity reflects changes in local neuronal circuit dysfunction and deficits in synaptic plasticity and 33 homeostasis that are widely propagated. 20 This is supported by our observation that the gene co-expression 34 module representing synaptic plasticity genes is downregulated cortex-wide and is significantly enriched with 35 common and rare ASD genetic risk variants, further emphasizing that synaptic plasticity is a convergent pathway 36 in ASD. Given this result, along with our observations of profound neuronal dysregulation present throughout the 37 ASD cortex, future work should determine which specific aspects of synaptic plasticity may contribute to causal 38 mechanisms in the disorder across specific brain regions and developmental timepoints. 39 Several additional factors should guide the interpretation of these results. The samples utilized in this 40 work were obtained from heterogeneous postmortem cortical tissue, meaning that the results reported here are 41 broadly applicable to the postnatal ASD cortex across both sexes and a span of ages from two to 68 years old, 42 and they should be interpreted in this context. Rigorous methodology was utilized at every step to account for 43 biological and technical variability, ensuring that the results reported here are conservative and widely applicable. 44 Additionally, bulk tissue RNA-seq, in contrast to single cell and nucleus RNA-seq, does not have the cellular 45 resolution to assess dissection variability across cortical regions and cell-type specificity of transcriptomic 46 changes. We addressed this by performing snRNA-seq, which significantly enhanced our understanding of 47 regional variation in ASD transcriptomic dysregulation. However, snRNA-seq also has its own limitations. While 48 snRNA-seq can profile tens of thousands of cells, snRNA-seq experiments typically have fewer unique samples 1 than bulk RNA-seq experiments, and the comparability of snRNA-seq cell-type proportions to true sample cell-2 type proportions is currently unclear. 23 It is also challenging to estimate isoform quantifications with single cell 3 RNA-seq approaches, whereas this remains a strength of bulk tissue RNA-seq. 24 Leveraging this, we 4 subsequently identified an upregulated isoform-specific co-expression module enriched with ASD GWAS 5 variants, implicating increased protein folding dysfunction for the first time as a putative pathway contributing to 6 ASD causal mechanisms. Interestingly, upregulated proteostasis is also implicated in Down's Syndrome, 25,26 7 supporting that protein folding machinery may be an affected biological process in multiple neurodevelopmental 8 disorders. The utilization of methods that have greater cellular resolution is necessary for the improved and 9 continued mapping of the results presented here to specific cortical cell-types. As we seek to gain a complete 10 understanding of ASD neural pathology, future approaches which integrate different sources of biological data -11 including this cortex-wide transcriptomic resource -to determine how ASD risk genes are acting in the brain will 12 be essential. 13 14 9

2
Sample Acquisition and Preparation for RNA-seq 3 4 Postmortem cortical brain samples were acquired from the Harvard Brain Bank as part of the Autism BrainNet 5 project (formerly the Autism Tissue Project, ATP) and the University of Maryland Brain Banks (UMDB Briefly, the limma-trend approach obtains normalized expression data through taking the log2(CPM) of read 1 counts with an adjustment for sample read depth variance. An offset value calculated with CQN 36 accounting 2 for GC content bias and gene/isoform effective length bias in read quantification was also incorporated during 3 the normalization process. With this normalized expression data, sample outliers were identified in each 4 sequencing batch by cortical lobe (frontal, parietal, temporal, and occipital) group that had both (1) an absolute 5 z-score greater than 3 for any of the top 10 expression principal components (PCs) and (2) a sample 6 connectivity score less than -2. Sample connectivity was calculated using the fundamentalNetworkConcepts 7 function in the WGCNA 10 R package, with the signed adjacency matrix (soft power of 2) of the sample biweight 8 midcorrelation as input. This process identified 34 outliers, resulting in a final total of 808 samples 9 (341=Control, 384=ASD, 83=dup15q) which were carried forward for analysis. 10 11 Evaluating Previous Co-Expression Modules and ASD DE Genes/Isoforms Cortex-wide 12 13 Linear models for all subsequent analyses are described in the Supplementary Methods. 14 15 To determine how gene co-expression modules previously identified in Parikshak et al. 5 and the residual. ASD-associated module eigengene region-specific ASD effects were identified using 21 contrasts (eg. Control_BA17 -ASD_BA17) with the limma 34 R package with this regressed expression dataset, 22 accounting for all biological covariates. Region-specific contrasts with a p-value < 0.05 were considered 23 significant (FDR-correction was unwarranted since only eight module eigengenes were examined). 24 25 To identify genes and isoforms dysregulated in ASD both within specific regions and cortex-wide, the limma 35 26 R package was applied with the gene and isoform expression data using our full gene and isoform models 27 (both biological and technical covariates). The standard limma 35 workflow was implemented as recommended 28 for linear mixed models. Region-specific dysregulation was identified as described above for the Parikshak et 29 al. 5 and Voineagu et al. 4 modules. Whole cortex dysregulation was established through subtracting the sum of 30 the ASD region-specific effects from the sum of the Control region-specific effects. For both region-specific and 31 whole cortex effects, genes and isoforms with an FDR-corrected p-value < 0.05 were considered significantly 32 dysregulated. dup15q region-specific and whole cortex dysregulation was also established in this manner. The 33 fixed effects of sex, age, and age 2 were also acquired (shared in Supplementary The methodology used to evaluate region-specific ASD effects compared to whole cortex ASD effects is 37 described in the Supplementary Methods. 38 39 Transcriptomic Regional Identity Analysis 40 41 To identify differentially expressed genes and isoforms between all 55 pairs of cortical regions, a regressed 42 gene expression dataset containing only the random effect of subject and the fixed effects of diagnosis and 43 region (along with the model residual) was used. Regression was performed as described for evaluation of 44 previously identified co-expression modules. Significant attenuation of DE genes between each pair of regions 45 (a reduction in transcriptomic regional identity differences) in ASD was established through the following 46 process. (1) ASD and Control subjects containing each region in the regional pair were extracted for use in the 47 analysis.
(2) Separately in ASD and Control subjects, the number of DE genes between regions was calculated 48 using the paired Wilcoxon signed-rank test. Genes with an FDR-corrected p-value < 0.05 were considered DE. 1 (3) The difference in the number of DE genes between regions for ASD v Control subjects was calculated (the 2 'true' difference). (4) A permuted distribution of the difference in DE genes between regions for ASD v Control 3 subjects was generated to test the 'true' difference. Each permutation (10,000 in total) randomly assigned 4 'ASD' and 'Control' status to subjects, but kept the number of ASD and Control subjects consistent with the 5 true number of ASD and Control subjects. (5) A two-tailed p-value was obtained from testing the 'true' 6 difference against the permuted distribution. If the regional comparison p-value < 0.05, with the number of DE 7 genes between regions in ASD less than that in Controls, then the regional comparison was considered 8 significantly attenuated in ASD. Otherwise, the regional comparison was considered over-patterned in ASD. 9 This procedure was repeated with isoform level regressed gene expression data (similarly, only containing the 10 random effect of subject and the fixed effects of diagnosis and region, along with the model residual) to identify 11 altered transcriptomic identities in ASD at the isoform-level. 12 13 The previously described permutation approach was designed to identify differences in transcriptomic regional 14 identity in ASD. Importantly, this method is not appropriate for assessing variance in expected numbers of DE 15 genes between regions across regional pairs and diagnoses, since the number of ASD and Control subjects 16 varied across regional pairs. To examine this, for each regional comparison we subset to 10 pairs of ASD and 17 Control subjects (10 was selected since every regional comparison had at least this many subjects). When 18 subsetting, subjects were removed such that the remaining subjects were closest in age to the median age of 19 the available samples for that regional comparison. A bootstrap approach was then used to calculate the 20 number of DE genes (p-value < 0.05) between regions separately in Control and ASD subjects through 21 sampling subjects with replacement (mean taken across 10,000 bootstraps). The same regressed expression 22 dataset used for the permutation approach was utilized for this bootstrap analysis. Any regional comparison in 23 which the number of DE genes between regions was less in ASD than in Control subjects was considered 24 trending towards attenuation in ASD. 25 26 To validate our bootstrapped estimates for the number of DE genes between pairs of regions in Controls, we 27 compared these estimates to those of the Allen Brain Atlas 9 , which is the best publicly available work for 28 comparison. Allen Brain Atlas regions were matched to Brodmann regions (Supplementary Table 4) and 29 matching regional pairs were extracted for comparison with this work. When the Allen Brain Atlas had two or 30 more regional pairs matching one regional pair in this work, the mean was taken across the Allen Brain Atlas 31 regional pairs. A p-value for the association of the number of DE genes between regions in Controls obtained 32 in this work compared to the Allen Brain Atlas was calculated from a linear model (cortex-wide bootstrap mean 33 ~ allen brain atlas mean). 34 35 We applied a stringent filtering process to identify high-confidence attenuated regional identity (ARI) genes 36 from each significantly attenuated regional comparison identified with the permutation procedure described 37 above. First, for each of the attenuated regional comparisons, we extracted the genes which were identified as 38 DE between regions in subjects labeled as Controls in each of the 10,000 permutations. Then, we calculated 39 how many times each of the genes truly DE between pairs of regions in the Control subjects were present in 40 their respective permuted groups (ranging from a possible 0 to 10,000 occurrences). Those 'true' DE genes 41 which were present in less than 95% of their respective permutations were retained as ARI genes for each 42 attenuated regional comparison. For each set of ARI genes (ten total), each gene was matched to the region in 43 which it had higher expression in Control subjects. The paired Wilcoxon signed-rank p-values identified for 44 these genes in Controls (those subjects used for the permutation analysis) were also extracted and are shared 45 in Supplementary Table 4. 46 47 ARI gene groups (ARI downregulated genes, those highly expressed in BA17 and BA39-40 relative to other 1 regions in Controls; ARI upregulated genes, those lowly expressed in BA17 and BA39-40 relative to other 2 regions in Controls) were created through taking the union (without duplicates) across all ten identified ASD-3 attenuated regional comparisons, and sorting genes into the two groups based on gene expression profiles 4 across regions. The details of this process are described in the Supplementary Methods, along with functional 5 annotation procedures. 6 7 Network-Based Functional Characterization 8 9 Standard workflows, as previously described in Parikshak et al. 5 and Gandal et al., 1 were followed (with minor 10 modifications) to identify gene and isoform co-expression modules using Weighted Gene Correlation Network 11 Analysis (WGCNA). 10 Details regarding network formation, module identification, and module functional 12 characterization are described in the Supplementary Methods. 13 14 snRNA-seq and Cell-type Deconvolution 15 16 Cell types were annotated based on expression of known marker genes visualized on the UMAP plot, violin 17 plots, and by performing unbiased gene marker analysis. To gain insight into the regional enrichment or 18 diagnostic enrichment of cell types, the relative proportion of the number of nuclei in each cell type was 19 normalized to the total number of nuclei captured from each library. Average cell-type proportions and standard 20 errors (across libraries) were scaled such that each Lobule x Diagnosis group sums to 100%, so that cell-type 21 proportions in these groups could be fairly compared across all cell-types. To determine if any changes in cell-22 type proportion were statistically significant, we implemented scDC 38 to bootstrap proportion estimates for our 23 samples (Supplementary Table 7). We employed a linear mixed model (random effect of subject) to 24 determine if any changes in cell-type proportion were present across regions and diagnoses. None of the 25 model covariates were statistically significant (p > 0.05 for all model covariates). However, we did find several 26 significantly different predicted cell-type proportions in ASD with cell-type deconvolution analysis. We describe 27 methods for cell-type deconvolution in detail in the Supplementary Methods. To identify genes differentially 28 expressed in ASD compared to control in each cell type, the non-parametric Wilcoxon rank sum test was 29 applied including gene detection rate and sequencing depth within the model. We compared frontal cortex ASD All of the raw bulk RNA-seq data (FASTQ files) and processed (utilized for DE gene analysis, transcriptomic 39 regional identity analysis, and WGCNA) bulk RNA-seq data that support the findings of this study will be 40 deposited in a publicly accessible repository. snRNA-seq data will be made available by the corresponding 41 authors upon reasonable request. All of the code, raw data, and processed data for the bulk RNA-seq analysis 42 that support the findings of this study will also be made available in a publicly accessible GitHub repository, 43 where readers may also access an R Shiny tool to visualize RNA-seq data across the biological covariates 44 assessed in this study. 45 46 47

Code availability
All of the code, raw data, and processed data for the bulk RNA-seq analysis that support the findings of this 3 study will also be made available in a publicly accessible GitHub repository, where readers may also access an 4 R Shiny tool to visualize RNA-seq data across the biological covariates assessed in this study. 5 6 b. Summary of sample composition (biological data, brain bank source, and PMI). testing the association of these covariates with module eigengenes is depicted. For the ASD and dup15q region-52 specific comparisons, cortical lobule colors are indicated (Fig. 1a), and bold-italic FDR p-values indicate that 53 these regions are effected with significantly greater magnitude than the ASD whole-cortex (Methods). For gene 54 biotypes, both positive and negative enrichment is shown (Methods). Positive enrichment is shown for cell-type  Top:  4 average-linkage hierarchical clustering of module eigengene biweight midcorrelations. Significant FDR corrected 5

Extended Data Figure Legends
p-values are indicated (FDR < 0.05; for GWAS, FDR < 0.1). Any signed -log10(p) colors greater or less than 5/-6 5 are set at a max/min of 5/-5 . For ASD, dup15q, and Age covariates, FDR p-value from the linear mixed model 7 testing the association of these covariates with module eigengenes is depicted. For the ASD and dup15q region-8 specific comparisons, cortical lobule colors are indicated (Fig. 1a). For gene biotypes, both positive and negative 9 enrichment is shown (Methods). Positive enrichment is shown for cell-type, GWAS, and rare variant enrichments 10 (Methods). Tissue, biological specimens or data used in this research were obtained from the Autism BrainNet (formerly 3 the Autism Tissue Program), which is sponsored by the Simons Foundation, and the University of Maryland 4 Brain and Tissue Bank, which is a component of the NIH NeuroBioBank. We are grateful to the patients and 5 families who participate in the tissue donation programs. Funding for this work was provided by grants to 6 D.H.G (NIMH R01MH110927, U01MH115746, P50-MH106438, and R01 MH-109912, R01 MH094714), grants 7 to MJG (SFARI Bridge to Independence Award, NIMH R01-MH121521, NIMH R01-MH123922, NICHD-P50-8 HD103557), grants to JRH (Achievement Rewards for College Scientists Foundation Los Angeles Founder 9 Chapter, UCLA Neuroscience Interdepartmental Program). We also thank Janet Sinsheimer for discussion of 10 the transcriptomic regional identity analysis methodology. 11 12  (8) 30 (5) 27 (5) 28 (8) 26 (5) 22 (7) 40 ( TFAP2A  HMOX1  OSR1  NR4A3  SNAI2  HES5  CREB3L1  TCF7  SOX11  ZBTB8B   0 10 20 (FDR) Methods overview for identifying statistically significant differences in transcriptomic regional identity in ASD. The regional comparison of BA17 v. BA41-42-22 is used here as an example. Number of DE genes between regions is calculated in controls and ASD samples (left). A permuted null distribution is then used to determine the significance of the difference in DE genes between controls and ASD samples (right). b. Regional comparisons with attenuation of transcriptomic regional identity in ASD with permutation p < 0.05 are connected with a bar (red scale), and a general trend towards cortex-wide attenuation is summarized by region color (blue scale; 0=no pairs exhibiting attenuation in ASD, 10=every pair exhibits attenuation in ASD). For region color, a regional pair is considered attenuated if it contains less DE in ASD compared to controls as measured with a bootstrap approach (Methods). Attenuated regional identity (ARI) genes are extracted from regional comparisons with permutation p < 0.05 (Methods). c-d. Overview of ARI downregulated (c.) and upregulated (d. test p-value is plotted above boxplots. Top right, PC 1 of ARI genes across all regions. Bottom left and bottom center, gene ontology and cell-type enrichment, respectively. Bottom right, top 10 attenuated transcription factors (TFs), where FDR is representative of how well these TFs distinguish BA17 and BA39-40 from the other nine cortical regions assessed here in controls (Methods). Enrichment for transcription factor binding sites is also depicted (Bonferroni-corrected p-value < 0.05 required for enrichment). Please see the Supplementary Methods for a description of all boxplots included in figures.