Creating and Benchmarking Spatial Dense Gene Expression Maps in the Human Cortex.

a, Spatially discontinuous Allen Human Brain Atlas (AHBA) microarray samples (red points) were aligned with MRI-derived cortical surface mesh reconstructions. b, AHBA vertex expression values were propagated using nearest-neighbor interpolation and subsequently smoothed (c). d, Subject-level maps were z-normalized and averaged to generate a single reference dense expression map (DEM) for each gene, as well as the associated expression gradient map (shown here for PVALB: top and bottom, respectively). e, DEMs can recover known expression boundaries in ISH data. Four canonical V1 area markers (Zeng et al., 2012 Cell) show a significantly sharp DEM expression gradient at the V1/V2 boundary (insert cortical map and Fig S2a,b), which is also evident in all four individual gene DEMs and DEM gradients (SYT6, PENK and Fig S2c). f, DEMs can discover previously unknown expression boundaries. Genes with high DEM gradients across the PeEc (parahippocampal) and TF (fusiform) gyri (inset cortical map) were validated in ISH data - showing sharp expression changes in both directions at this boundary (CHRNA3, NGB and Fig Sd-f).g, Illustrative comparisons of selected DEMs against regional variation in microscale measures of cellular composition : scatterplot showing the global correlation of regional cellular proportions from single nucleus RNAseq (snRNAseq) across 16 cells and 6 regions(Lake et al., 2016) with DEM values for corresponding cell-type marker genes (R=0.48, pspin<0.001, excluding Ex3-V1 and In8-BA10 outlier samples). h, DEMs for markers of 6 neuronal subtypes (3 excitatory: FEZF2, RORB, THEMIS, 3 inhibitory: PVALB, SST, VIP) based on recently validated subtype marker genes(Bakken et al., 2021; Hodge et al., 2019)i, Illustrative comparison of layer IV marker DEMs with corresponding mesoscale cortical measure of layer IV thickness from a 20μm 3D histological atlas of cortical layers. j, Illustrative comparisons of selected DEMs with corresponding macroscale cortical measures from independent neuroimaging markers.

Mapping transcriptional distinctiveness in the human cortex and its alignment with macroscale structure and function.

a, Regional transcriptomic distinctiveness (TD) can be quantified as the mean absolute z-score of dense expression map (DEM) values at each vertex (top), and visualized as a continuous cortical map (middle, TD encoded by color) or in a relief map of the flattened cortical sheet (bottom, TD encoded by color and elevation, Sup Movie 1). Black lines on the inflated view identify cuts for the flattening procedure. The cortical relief map is annotated to show the central sulcus (CS), and peaks of TD overlying dorsal sensory and motor cortices (Brodmann Areas, BA2, BA4), the primary visual cortex (V1), temporal pole (TGd), insula (Ins) and ventromedial prefrontal cortex (OFC). b, Thresholding the TD map through spatial permutation of DEMs (tspin Methods) and clustering significant vertices by their expression profile defined six TD peaks in the adult human cortex (depicted as coloured regions on terrain and inflated cortical surfaces). c, Cortical vertices projected into a 3D coordinate system defined by the first 3 principal components (PCs) of gene expression, coloured by the continuous TD metric (left) and TD peaks (right). TD peaks are focal anchors of cortex-wide expression PCs d, TD peaks show statistically-significant functional specializations in a meta-analysis of in vivo functional MRI data. e, The average magnitude of local expression transitions across genes (color) and principal orientation of these transitions (white bars) varies across the cortex. f, Cortical folds in AHBA donors (top surface maps and middle flat-map) tend to be aligned with the principal orientation of TD change across cortical vertices (p<0.01, middle histogram, sulci running perpendicular to TD change), and the strength of this alignment varies between cortical regions. g, Putative cortical areas defined by a multimodal in vivo MRI parcellation of the human cortex(Glasser et al., 2016) (top surface maps and middle flat-map) also tend to be aligned with the principal direction of gene expression change across cortical vertices (p<0.01, middle histogram, sulci running perpendicular to long axis of area boundaries), and the strength of this alignment varies between cortical areas.

Cortex-wide Gene Coexpression Patterns Reflect Multiple Spatial Scales and Developmental Epochs of Brain Organization.

a, Overview of Weighted Gene Co-expression Network Analysis (WGCNA) pipeline applied to the full DEM dataset. Starting top left: the pairwise DEM spatial correlation matrix is used to generate a topological overlap matrix between genes (middle top) which is then clustered. Of the 23 WGCNA-defined modules, 7 were significantly enriched for non-cortical genes and removed, leaving 16 modules. Each module is defined by a set of spatially co-expressed genes, for which the principal component of expression can be computed and mapped at each cortical point (eigenmap). M6 is shown as an example projected onto an inflated left hemisphere (M6 z-scored expression and M6 expression change), and the bulk transcriptional distinctiveness (TD) terrain view from Fig 2 (M6 expression). b, The extremes of WGCNA eigenmaps highlight different peaks in the cortical terrain: the main TD terrain colored by TD value (center, from Fig 2), surrounded by TD terrain projections of selected WGCNA eigenmaps. c, WGCNA modules (eigenmaps and gradient maps, rows) are enriched for multiscale aspects of cortical organization (columns). Cell color intensity indicates pairwise statistical significance (p<0.05), while black outlines show significance after correction for multiple comparisons across modules. Columns capture key levels of cortical organization at different spatial scales (arranged from macro- to microscale) and developmental epochs: spatial alignment between module eigenmaps and in vivo MRI maps of cortical folding orientation, cortical thickness and T1/T2 ratio, fMRI resting-state functional networks; enrichment for module gene sets for independent annotations (Table S2) marking: cortical layers(He et al., 2017; Maynard et al., 2021); cell types(Darmanis et al., 2015; Habib et al., 2017; Hodge et al., 2019; Lake et al., 2018, 2016; Li et al., 2018; Ruzicka et al., 2021; Velmeshev et al., 2019; Zhang et al., 2016); subcellular compartments(Binder et al., 2014); synapse-related genes(Koopmans et al., 2019); protein-protein interactions between gene products (Szklarczyk et al., 2019); temporal epochs of peak expression(Werling et al., 2020) [“fetal”: 8-24 21 post conception weeks (PCW) / “perinatal’’ 24 PCW-6 months / “postnatal” >6 months]; transient layers of the mid-fetal human cortex at 21 post conception weeks (PCW)(Miller et al., 2014)[subpial granular zone (SG), marginal zone (MZ), cortical plate (CP), subplate (SP), intermediate zone (IZ), subventricular zone (SZ) and ventricular zone (VZ)]; and fetal cell types at 17-18 PCW(Polioudakis et al., 2019). d, Independent validation of multiscale enrichments for selected modules M2 & M12. M2 significantly overlaps the Neurosynth topic associated with the terms motor, cortex and hand. Two high-ranking M2 genes, MOG & TF exhibit clear layer VI peaks on ISH and GO enrichment analysis myelin-related annotations. M12, overlapping the limbic network most closely overlapped the Neurosynth topic associated with social reasoning. Two high-ranking M22 genes GABRA2 and GRIN2B showed layer II ISH peaks and GO enrichment analysis revealed synaptic annotations. e, Network visualization of pairwise overlaps between annotational gene sets used in Fig 3c, including WGCNA module gene sets (inset expression eigenmaps).

ASD risk genes follow two different spatial patterns of cortical gene expression which differentially predict cortical changes in ASD.

a, Enrichment of WGCNA module gene sets for risk genes associated with atypical brain development through enrichment of rare deleterious variants in studies of Autism Spectrum Disorder (ASD), Schizophrenia (Scz), severe developmental disorders (DDD, Deciphering Developmental Disorders study) and Epilepsy. Cell color intensity indicates pairwise statistical significance (p<0.05), while outlined matrix cells survived correction for multiple comparisons across modules. b, Summary of multiscale and developmental annotations from Fig 3c for M12 and M15: the only two WGCNA modules enriched for risk genes of more than one neurodevelopmental disorder. c, M12 and M15 genes clustered by the strength of their membership to each module. Color encodes module membership. Shape encodes annotations for two GO Biological Process annotations that differ between the module gene sets: neuronal communication and regulation of gene expression. Text denotes specific ASD risk genes. d, contrasting GO enrichment of M12 and M15 for neuronal communication and regulation of gene expression GO Biological Process annotations. e, M12 and M15 differ in the developmental trajectory of their average cortical expression between early fetal and mid-adult life(Li et al., 2018). f, Regional differences in intrinsic expression of the M15 module (but not the M12 module) in adult cortex is correlated with regional variation in the severity of altered cortical gene expression (number of differentially expressed genes) in ASD(Haney et al., 2020). g, Statistically-significant regional alterations of cortical thickness (CT) in ASD compared to typically developing controls from in vivo neuroimaging(Di Martino et al., 2017, 2013) (top). Areas of cortical thickening show a statistically-significant spatial overlap (Dice overlap = 0.68, pspin<0.01) with regions of peak intrinsic expression for M15 in adult cortex (bottom). h, M15 eigenmap expression (but not M12 eigenmap) shows significant spatial correlation with relative cortical thickness change in ASD.

Statistical tests used to compare spatial maps and gene sets derived from the Allen Human Brain Atlas with independent multiscale neuroscientific resources.