Prioritizing Genetic Contributors to Cortical Alterations in 22q11.2 Deletion Syndrome Using Imaging Transcriptomics

Abstract 22q11.2 deletion syndrome (22q11DS) results from a hemizygous deletion that typically spans 46 protein-coding genes and is associated with widespread alterations in brain morphology. The specific genetic mechanisms underlying these alterations remain unclear. In the 22q11.2 ENIGMA Working Group, we characterized cortical alterations in individuals with 22q11DS (n = 232) versus healthy individuals (n = 290) and conducted spatial convergence analyses using gene expression data from the Allen Human Brain Atlas to prioritize individual genes that may contribute to altered surface area (SA) and cortical thickness (CT) in 22q11DS. Total SA was reduced in 22q11DS (Z-score deviance = −1.04), with prominent reductions in midline posterior and lateral association regions. Mean CT was thicker in 22q11DS (Z-score deviance = +0.64), with focal thinning in a subset of regions. Regional expression of DGCR8 was robustly associated with regional severity of SA deviance in 22q11DS; AIFM3 was also associated with SA deviance. Conversely, P2RX6 was associated with CT deviance. Exploratory analysis of gene targets of microRNAs previously identified as down-regulated due to DGCR8 deficiency suggested that DGCR8 haploinsufficiency may contribute to altered corticogenesis in 22q11DS by disrupting cell cycle modulation. These findings demonstrate the utility of combining neuroanatomic and transcriptomic datasets to derive molecular insights into complex, multigene copy number variants.


Introduction
22q11.2 deletion syndrome (22q11DS) arises from the deletion of a segment of chromosome 22 due to misalignment of low copy repeats (LCR) during nonallelic homologous recombination. It occurs in approximately 1 in 3000-4000 births and spans a ∼2.6 megabase (Mb) region that results in the hemizygous deletion of 46 protein-coding genes in 85-90% of patients (Guna et al. 2015), with ∼10-15% of 22q11DS patients carrying a smaller, nested deletion (McDonald-McGinn et al. 2015). 22q11DS is associated with a broad phenotype that includes heart anomalies, immune dysfunction, and high rates of neuropsychiatric and neurodevelopmental disorders such as schizophrenia, intellectual disability, and autism spectrum disorder (ASD; Jonas et al. 2014;Schneider et al. 2014). Alterations in brain structure and function are thought to contribute to the psychiatric and developmental phenotypes frequently observed in the disorder.
Indeed, it is now established that 22q11DS is associated with widespread alterations in brain morphology. Early magnetic resonance imaging (MRI) studies reported whole brain volumetric reductions in 22q11DS, with greater reductions in midline regions, as well as in posterior relative to anterior regions (Tan et al. 2009;Karayiorgou et al. 2010). However, as cortical gray matter volume reflects the product of cortical surface area (SA; i.e., area covered by the cortex) and cortical thickness (CT; i.e., thickness of the 6 neocortical layers), which appear to be determined through relatively independent genetic and neurodevelopmental mechanisms (Panizzon et al. 2009;Winkler et al. 2010;Chen et al. 2013;Grasby et al. 2020, but see also Schmitt et al. 2018), recent studies have examined these morphometric characteristics separately. Thus, brain volume reductions in 22q11DS were recently found to be driven by widespread reductions in SA (Sun et al. 2019). Conversely, CT tends to be increased in 22q11DS, with focal thinning in only a minority of regions. Overall, the magnitude of SA alterations in 22q11DS is roughly 2-fold the magnitude of CT alterations . Importantly, while 22q11.2 deletions yield reduced expression of the majority of genes within the locus (Stark et al. 2008;Jalbrzikowski et al. 2015;Lin et al. 2016;Gordon et al. 2019), all genes within the locus are not expected to contribute equally to brain phenotypes in the disorder (Motahari et al. 2019). As neuroanatomic abnormalities are associated with a range of neuropsychiatric and developmental phenotypes, clarifying the individual genes underlying these abnormalities may provide insight into molecular mechanisms that contribute to broader psychiatric and developmental phenotypes in 22q11DS.
Leveraging comprehensive maps of gene expression in the human brain offers one promising approach to identify molecular mechanisms underlying neuroanatomic deviations in 22q11DS. Recent studies have used the Allen Human Brain Atlas (AHBA), a transcriptomic dataset quantifying the expression of over 20 000 genes across postmortem brain tissue from six psychiatrically healthy individuals, to elucidate mechanisms underlying cellular and neural circuit variation in healthy individuals and in populations with neuropsychiatric and neurodegenerative disorders (Fornito et al. 2019). By examining the spatial convergence between brain phenotypes and gene expression patterns, recent studies found that brain regions that are closer in physical proximity (Hawrylycz et al. 2012) or have functionally correlated activity (Richiardi et al. 2015) share more similar transcriptomic expression patterns. Similarly, a prominent rostro-caudal gradient of gene expression has been found across the cortex (Bernard et al. 2012;Hawrylycz et al. 2012;Miller et al. 2014), which is thought to reflect the rostrocaudal gradient of neurogenesis and cell composition in which posterior brain regions have a higher density of neurons that are smaller in size, while anterior regions tend to have a lower density of neurons that are larger in size and spine density (Cahalane et al. 2012;Charvet et al. 2015;Fornito et al. 2019). In clinical populations, the pattern of structural dysconnectivity in schizophrenia patients was spatially correlated with the expression of 43 genes previously implicated in schizophrenia by genome-wide association (Romme et al. 2017); regional expression of the Parkinson's risk gene, MAPT, was spatially correlated with the topography of connectivity differences in patients with Parkinson's disease ); and regional expression of transcriptionally down-regulated genes in postmortem cortex of ASD patients was associated with severity of CT deviation in ASD (Romero-Garcia 2019). Spatial convergence analyses were also recently applied in a 16p11.2 deletion mouse model to identify genes within the locus that may be causally related to structural brain changes associated with the copy number variant (CNV; Kumar et al. 2018). Thus, prior studies of neurodevelopmental, neuropsychiatric, and neurodegenerative populations suggest that identifying genes with expression patterns that are spatially correlated with neuroimaging phenotypes can offer a useful strategy to elucidate genetic drivers of altered brain structure and function.
Here, as part of the 22q11.2 Enhancing Neuroimaging Genetics through Meta-Analysis (ENIGMA) Working Group (Thompson et al. 2020), we therefore integrated neuroanatomic data from a large multicenter cohort of 22q11DS individuals with molecularly confirmed deletions spanning the full LCR A-D region and transcriptomic data from the AHBA. By characterizing the spatial convergence between regional expression of each individual 22q11.2 gene and the severity of morphometric alterations within patients, we sought to systematically prioritize individual genes within the 22q11.2 locus that may be causally related to these alterations and elucidate potential underlying molecular mechanisms.

Structural MRI Data
Structural MRI (sMRI) data from 386 22q11DS patients and 315 typical developing controls analyzed in a previously published study from the 22q11DS ENIGMA working group  were used to derive measures of SA and CT deviance in 22q11DS for the current study. Briefly, in the original study, data were pooled across nine study sites with patient and control data. FreeSurfer image processing software (version 5.3.0; http://surfer.nmr.mgh.harvard.edu) was used to process 1 mm 3 T 1 -weighted structural images acquired with an MPRAGE sequence. Quality control was implemented using validated and standardized processing pipelines developed for the ENIGMA consortium (Thompson et al. 2014(Thompson et al. , 2017http://enigma.ini.u sc.edu/protocols/imaging-protocols). Total intracranial volume (ICV) and SA and CT measures for 68 cortical regions (34 per hemisphere) were calculated based on the Desikan-Killiany atlas. Group effects in this multisite 22q11DS cohort were previously found to be highly consistent across sites .
The purpose of the current analysis was to examine whether the regional expression patterns of individual genes in the 22q11.2 locus are associated with regional severity of SA or CT deviance in 22q11DS patients. Consequently, our analyses focused on a homogeneous sample of 22q11DS patients with the full ∼2.6 Mb A-D deletion and the typical expression patterns (i.e., in healthy individuals) of the corresponding genes within this region, based on the AHBA. 22q11.2 deletion breakpoints for each patient were determined using multiplex ligationdependent probe amplification (MLPA; Vorstman et al. 2006), which is a polymerase chain reaction (PCR)-based assay that can detect copy number deletions and duplications for up to 50 DNA probe sequences in one reaction. Due to its low cost, high sensitivity and specificity, and medium throughput, it is considered a gold standard method for CNV genotyping in humans (Kerkhof et al. 2017). MLPA for the current study was completed using the SALSA MLPA Probemix P250-B2 DiGeorge kit from MRC-Holland, which includes 29 probe sequences within the 22q11.2 locus to discern between common 22q11.2 deletion subtypes. Sites or scanners with no 22q11DS patients with a confirmed A-D deletion or corresponding healthy control data were excluded from the current analyses, leaving data from 232 22q11DS patients with confirmed A-D deletions and 290 controls collected on 10 scanners across 8 sites for analysis. Regional SA measures were adjusted for effects of age, sex, and site/scanner; different scanners were treated as independent "sites." Regional CT measures were additionally adjusted for age 2 , based on significant nonlinear effects of age previously found in most ROIs for CT .
Subject consent at each site was obtained according to the Declaration of Helsinki, and study protocols were approved by ethical committees at each institution. Detailed information on the sample recruitment procedures, image acquisition parameters, and data processing are published elsewhere ).

Gene Expression Data
AHBA transcriptomic data for 20 737 largely protein-coding genes, registered to the Desikan-Killiany cortical atlas for integration with FreeSurfer-based analyses, were obtained from https://figshare.com/articles/A_FreeSurfer_view_of_the_cortica l_transcriptome_generated_from_the_Allen_Human_Brain_Atla s/1439749 (French and Paus 2015). The original AHBA assayed the expression of 58 692 probes at a high spatial resolution, using custom Agilent arrays in 3702 brain samples derived from six healthy adults with no known neuropsychiatric or neuropathological history (Hawrylycz et al. 2012). In the French and Paus (2015) atlas, expression values from multiple probes for a given gene were first averaged to yield one expression value per gene per tissue sample. Each cortical brain tissue sample was then mapped to the nearest Desikan-Killiany cortical region based on its Montreal Neurological Institute coordinates, and for each individual brain, the median expression value across tissue samples mapping to a given Desikan-Killiany region was calculated for each gene. The median expression level across the six brains, per region, was then calculated for each of the 20 737 genes (see French and Paus 2015 for details). Because expression levels expression levels were measured from left hemisphere regions in all six brains, but only in two brains for right hemisphere regions, all region-based analyses used only left hemisphere regions.
To define protein-coding genes within the 22q11.2 locus, coordinates for the 22q11.2 locus were obtained from genomewide studies of CNVs associated with schizophrenia (Marshall et al. 2017) and ASD (Sanders et al. 2015). Marshall et al. (2017) reported CNV borders in hg18; the UCSC LiftOver tool was used to convert them to hg19. As the 22q11.2 locus defined in these studies shared more than 90% overlap in basepairs (bp), the final 22q11.2 boundaries were defined as the union between those identified in these two studies . HGNC gene symbols for protein-coding genes within the locus were retrieved from Ensembl using the BioMart package in R (Durinck et al. 2009). Genes with mean log 2 expression levels > 5, averaged across brain regions, were considered brain expressed, leaving 28 brain-expressed 22q11.2 protein-coding genes for investigation. The consistency of the expression of each 22q11.2 gene across the six donors, as defined by French and Paus (2015), is reproduced in the Supplementary Information.

Demographic Characteristics
Group differences in age and sex were assessed with a univariate ANOVA and a chi-squared test, respectively.

Neuroanatomic Group Differences
To characterize neuroanatomic alterations in the 22q11DS patients included in the present analyses, global and regional SA and CT metrics were first adjusted for age, sex, and scanner effects by conducting linear models with each neuroanatomic metric set as the dependent variable and with age, sex, and scanner set as independent variables, additionally including age 2 for CT measures . Residualized values for each subject and neuroanatomic metric were retained from the linear models. Group differences in covariate-adjusted total SA and in mean CT were then examined using general linear models with total SA or mean CT as the dependent variable and group as the independent variable. Group differences in covariate-adjusted SA and CT for each left hemisphere region were similarly tested using general linear models. This analysis flow was selected such that between-group neuroanatomic analyses and spatial convergence analyses were conducted on the same covariate-adjusted neuroanatomic values. Group differences in regional neuroanatomic measures are shown corrected for multiple comparisons using false discovery rate (FDR) correction (q-value < 0.05 across 34 regions).

Prioritizing 22q11.2 Genes Based on Spatial Convergence of Gene Expression and Neuroanatomic Deviance
To prioritize 22q11.2 genes that may be causally involved in SA or CT alterations in 22q11DS, our primary analysis focused on correlations between regional expression of each brain-expressed, protein-coding 22q11.2 gene and regional deviance in SA and CT, respectively. The neuroanatomic deviance of 22q11DS patients compared with controls was first defined using Z-scores for each covariate-adjusted regional measure of SA and CT. Normalized deviance scores were utilized over raw group difference scores in order to account for differences in the area or average thickness of regions as they are defined in the Desikan-Killiany atlas (i.e., to avoid nonmeaningful larger deviance scores in regions that are defined in the reference Desikan-Killiany atlas as covering a greater number of voxels; Grothe et al. 2018). Thus, mean SA and CT for each left hemisphere region were first calculated for each group. Given the global tendency for 22q11DS patients to show smaller SA overall, mean SA per region for 22q11DS patients was subtracted from the mean for controls and divided by the standard deviation for controls for each region to yield a Z-score severity measure of 22q11DS SA deviance ( SA) per region. For CT, given the global tendency for 22q11DS patients to show higher CT, mean CT per region for controls was subtracted from mean CT scores for 22q11DS patients and divided by the standard deviation for controls per region to yield a Z-score severity measure of 22q11DS CT deviance ( CT) per region.
Pearson's correlation coefficients were then used to examine spatial convergence in the expression of each brain-expressed, protein-coding AHBA gene and regional variation in 22q11DS SA and CT severity. Pearson's correlations were used given that SA, CT, and the expression of the majority of brainexpressed, protein-coding AHBA genes (>75%) were normally distributed across regions (Shapiro-Wilk Test P > 0.05) and that these are continuous variables. Spearman's (nonparametric) correlations yielded highly similar results (see Supplementary  Information). Error in correlation coefficient estimates and 95% confidence intervals were assessed with bootstrapping (i.e., resampling the 34 cortical regions 1000 times with replacement) using the "boot" package in R. The ratio of the mean correlation per gene to its bootstrap standard deviation (i.e., Z-score correlation) was used to generate percentile ranks for all brainexpressed, protein-coding genes indexed in AHBA (i.e., 10 344 genes) and to derive corresponding P-values based on this empirical distribution of all Z-score correlations across the AHBA. 22q11.2 genes with expression patterns that showed positive spatial correlations with 22q11DS SA or CT severity at extreme high-rank values relative to all brain-expressed, protein-coding AHBA genes (P AHBA < 0.05) were considered statistically significant Seidlitz et al. 2020).
This analysis leverages the fact that the proximal consequence of the 22q11.2 deletion is reduced expression of the majority of deleted genes and assumes that greater severity of neuroanatomic deviation in 22q11DS will be evident in regions where potential causally related genes within the locus are typically most highly expressed.

Prioritizing 22q11.2 Genes Based on Partial Least Squares Regression
As a complementary approach to the above described primary analysis that leveraged each individual gene and individual SA and CT deviance score per region, we also examined whether an alternate approach utilizing a partial least squares regression (PLSR) data reduction technique would prioritize similar 22q11.2 genes. This secondary analysis was implemented to test the robustness of our primary findings to variation in analysis techniques. Thus, PLSR identifies principal components based on both the predictors (i.e., here, the expression of all brainexpressed, protein-coding genes in AHBA) and the outcome (i.e., here, 22q11DS SA or CT severity) to maximally explain covariance between predictors and the outcome (Wehrens and Mevik 2007). In imaging transcriptomic analyses utilizing PLSR, the first principal component (PLS1) represents the linear combination of gene weights with expression patterns that best predict the neuroimaging measure across regions. PLSR has previously been used to identify biological processes that are broadly associated with spatial neuroanatomic deviance patterns in neuropsychiatric and neurodevelopmental populations (Whitaker et al. 2016;Seidlitz et al. 2018;Morgan et al. 2019) and is generally used in scenarios when prioritizing individual genes out of an a priori set of candidate genes is not the primary analysis goal. Nevertheless, gene weights on the first principal component of a given model can also be used to rank genes with expression patterns that best predict neuroanatomic alterations.
To examine the validity of PLS1 for our SA and CT models, the significance of the variance explained by PLS1 for each model was tested by permuting the outcome labels 10 000 times. To establish gene rankings for PLS1 for each model, error in estimating each gene's PLS1 weight for the SA and CT models was assessed by bootstrapping (i.e., resampling the 34 cortical regions 1000 times with replacement). The ratio of the mean loading weight of each gene to its bootstrap standard deviation (i.e., Z-score loading) was used to generate percentile ranks for all genes for their PLS1 loading relative to the distribution across all brain-expressed, protein-coding genes in AHBA and derive corresponding empirical P-values. Genes with extreme highrank values (P AHBA < 0.05) within this empirical distribution were considered to load significantly on PLS1 for each model. This empirical distribution was also used to prioritize 22q11.2 genes within the PLSR approach. Finally, gene ontology (GO) analyses examined whether all genes that loaded significantly on PLS1 for each model were enriched for specific biological pathways, molecular functions, or cellular components using g:Profiler (Reimand et al. 2016), with "moderate" hierarchical filtering (best per parent) and a minimum query/term overlap size of 5 genes. Only pathways with 10 to 2000 genes were included, and a custom background was set to all protein-coding, brain-expressed AHBA genes.

Characterizing 22q11.2 Gene Prioritization Relative to Top Genes in Random Gene-Lists
Given that the proximal consequence of 22q11.2 deletions is reduced expression of genes within the locus, our analyses focused on prioritizing individual genes within the locus that are most likely to be causally related neuroanatomic alterations in 22q11DS. Nevertheless, to contextualize the strength of the correlations between the top 22q11.2 genes and SA or CT severity in 22q11DS, we also generated 10 000 random lists of 28 brain-expressed, protein-coding AHBA genes. The correlation Z-scores for the top 22q11.2 gene identified in the primary 22q11SDS SA or CT severity analyses, as well as the Z-score loadings for the top 22q11.2 gene identified in the PLSR analyses  (198219) for SA or CT severity, were compared with the distribution of top correlation Z-scores or top PLSR Z-score loadings across the 10 000 lists of 28 random genes (P GENE-LIST ).

Deficiency-Induced Down-Regulated miRNAs
Given the prominent gene regulatory role of DGCR8 via microRNA (miRNA) biogenesis, to follow up on spatial convergence results between regional expression of DGCR8 and 22q11DS cortical SA severity, we characterized the gene targets of miRNAs previously suggested to be down-regulated in the cortex due to DGCR8 deficiency. Specifically, given that cortical tissue derived directly from 22q11DS patients was not available, we focused on miRNAs with significantly reduced expression in prefrontal cortex (PFC) in a mouse model of 22q11DS, whose down-regulation was accounted for by DGCR8 deficiency (Stark et al. 2008). Down-regulated miRNA names were converted from miRBase version 9.1 to miRBase version 21.0 nomenclature using miRNA Accession IDs, and the human gene targets of the homologous human miRNAs were identified using miRTarBase v7.0 (Chou et al. 2018). Gene targets of DGCR8 deficiency-induced down-regulated miRNAs were functionally annotated using GO biological pathways, molecular functions, and cellular components from g:Profiler (Reimand et al. 2016). Gene targets were also tested for enrichment for lists of genes expressed in specific cell types and specific human brain regions during specific developmental periods (i.e., relative to all other regions/developmental periods) using the Specific Expression Analysis tool (http://genetics.wustl.edu/jdlab/cseatool-2/; Dougherty et al. 2010). See Supplementary Material for details.
In line with the larger 22q11.2 ENIGMA study  SA reductions in 22q11DS were widespread, with particularly prominent reductions in midline posterior brain regions, including the cuneus, precuneus, and lingual gyrus, as well as lateral association regions including superior parietal cortex and rostral middle frontal gyrus ( Fig. 1A; Supplementary Table 1). Parallel analyses adjusting for ICV identified similar regions of lower SA in 22q11DS versus controls (Supplementary Table 2).
Regional CT differences were also similar to those found in the previous 22q11.2 ENIGMA study, with the majority of regions showing subtle increases in CT in 22q11DS that were greatest in frontal and parietal regions including rostral and caudal middle frontal gyrus, medial and lateral orbitofrontal cortex, precentral and postcentral gyrus, and supramarginal gyrus, as well as in pericalcarine cortex and insula ( Fig. 1D; Supplementary  Table 3). Significant focal thinning in 22q11DS was found in the caudal anterior cingulate, superior temporal cortex, and parahippocampus.

Prioritized Genes Associated with Neuroanatomic Alterations in 22q11DS
Spatial correlation analyses comparing the expression patterns of 22q11.2 genes to SA severity in 22q11DS revealed a significant positive correlation between regional expression of DGCR8 and regional SA severity, Pearson r = 0.53, and P AHBA = 0.006 ( Fig. 1A-C; Table 2). AIFM3 expression was additionally associated with SA severity, Pearson r = 0.42, and P AHBA = 0.041. Thus, brain regions with higher expression of DGCR8 and AIFM3 in healthy individuals showed greater reductions in cortical SA in 22q11DS (Supplementary Fig. 1). Parallel analyses using SA measures additionally adjusted for ICV, as well as those using Spearman's correlations, yielded highly similar results (Supplementary Tables 4 and 5, respectively).
Parallel analyses for CT revealed that the spatial pattern of CT severity in 22q11DS was significantly associated with regional expression of P2RX6, r = 0.43, and P AHBA = 0.022 ( Fig. 1D-E; Table 3; Supplementary Fig. 2). Thus, regions with Figure 1. Variation across 34 left hemisphere cortical regions in: A) Z-score surface area deviance ( SA) severity in 22q11DS patients relative to controls (higher Z-score indicates region with greater reduction in SA in 22q11DS); B) DGCR8 expression; C) AIFM3 expression; D) Z-score cortical thickness deviance ( CT) severity in 22q11DS relative to controls (higher Z-score indicates region with greater increase in CT in 22q11DS); and E) P2RX6 expression. Expression of DGCR8 and AIFM3 were significantly associated with SA in 22q11DS and expression of P2RX6 was significantly associated with CT in 22q11DS.
higher P2RX6 expression in healthy individuals showed greater increases in CT among 22q11DS patients compared with healthy controls. Analyses using Spearman's correlations yielded similar results (Supplementary Table 6).
Restricting analyses to only brain-expressed, protein-coding genes that were consistently expressed across the 6 AHBA donors (average donor-to-median expression ρ > 0.446, as defined by French and Paus (2015); 4947 genes) included 16 22q11.2 genes and similarly identified significant associations between regional expression of DGCR8 and SA severity in 22q11DS (Supplementary Table 7) and regional expression of P2RX6 with CT severity (Supplementary Table 8).

Robustness of Neuroanatomic Deviance in 22q11DS and Prioritized 22q11.2 Genes Across Age Subgroups
To examine the robustness of these relationships across development, we carried out parallel analyses examining group differences and gene expression spatial correlations with SA and CT severity for 3 age subgroups in the 22q11.2 ENIGMA cohort: children (≤12 years; 22q11DS n = 75, HC n = 89), adolescents (13-17 years; 22q11DS n = 67, HC n = 77), and adults (≥18 years; 22q11DS n = 90, HC n = 124). The overall pattern of neuroanatomic differences among 22q11DS patients was highly similar for each age subgroup compared with the overall 22q11DS sample (see Supplementary Tables 9 and 10, respectively). Consistent with this, conducting an exploratory omnibus linear model across age groups, including a CNV group by age group interaction term revealed no significant CNV group by age group interactions for any ROI (FDR corrected P-values > 0.05), similar to findings in the larger ENIGMA cohort ). This suggests that neuroanatomic differences in SA and CT in 22q11DS patients are largely established by childhood.

Similar Prioritized 22q11.2 Genes Identified Using Partial Least Squares Regression
Prioritizing 22q11.2 genes based on gene loadings on PLS1 for the 22q11DS SA and CT models highlighted similar genes. Thus, PLS1 for the SA model explained 26.7% of the covariance between regional 22q11DS SA severity and gene expression, which was significantly more than expected by chance (permutation P = 0.002). Among the 22q11.2 genes, DGCR8 and AIFM3 loaded significantly on PLS1 (P AHBA < 0.05; Supplementary Table 13). Functional annotation of all genes with significant loadings on PLS1 indicated that the spatial pattern of SA in 22q11DS was broadly associated with gene regulatory processes, including transcriptional activity and chromatin organization ( Supplementary Fig. 3A).
PLS1 for the CT model explained 23.9% of the covariance between regional 22q11DS CT severity and gene expression, which was more than expected by chance (permutation P = 0.003). Among 22q11.2 genes, P2RX6 loaded significantly on PLS1 (P AHBA = 0.007; Supplementary Table 14). Functional annotation of all genes with significant loadings on PLS1 indicated that the spatial pattern of CT in 22q11DS was broadly associated with genes involved in transmembrane and ion transmembrane transport ( Supplementary Fig. 3B).

Characterizing 22q11.2 Gene Prioritization Relative to Random Genome-Wide Gene-Lists
Relative to 10 000 lists of 28 random brain-expressed, proteincoding AHBA genes, the strength of the relationship between 22q11DS SA severity and the correlation Z-score for the top 22q11.2 gene, DGCR8, was greater than the strongest association observed in 85.9% of random gene-lists, P GENE-LIST = 0.141. For 22q11DS CT severity, the strength of the relationship with the top 22q11.2 gene, P2RX6, was greater than 63.4% of lists, P GENE-LIST = 0.366. The top 22q11.2 gene loading on PLS1 for the 22q11DS SA model, DGCR8, was similarly greater than for 94.8% of randomly generated gene-lists, P GENE-LIST = 0.052, and the top 22q11.2 gene loading for PLS1 for the 22q11DS CT model, P2RX6, was greater than for 83.4% of gene-lists, P GENE-LIST = 0.166. The elevated but still modest ranking for the strongest 22q11.2 gene spatial correlation with 22q11DS SA and CT severity, respectively, compared with top correlations identified among random sets of any 28 AHBA genes likely reflects the fact that groups of genes work in concert to carry out specific biological functions and, relatedly, that the spatial expression patterns of many genes are highly correlated (Hawrylycz et al. 2012).

Characterizing Downstream Consequences of DGCR8 Haploinsufficiency
DGCR8 was the most prominent gene within the 22q11.2 locus associated with SA severity in 22q11DS. Notably, DGCR8 is a core component of the miRNA microprocessor complex involved in the biogenesis of miRNAs, which are small noncoding RNAs that critically regulate gene expression by binding target messenger RNA (mRNA) transcripts to accelerate their degradation or silence their translation into proteins. Given this prominent gene regulatory role of DGCR8, to better understand mechanisms through which DGCR8 haploinsufficiency may contribute to SA in 22q11DS, we therefore explored whether gene targets of miRNAs that are down-regulated in mouse PFC due to DGCR8 deficiency (Stark et al. 2008) converge on specific biological processes. The 59 miRNAs down-regulated in PFC (Supplementary Table  15) due to DGCR8 deficiency targeted 6804 unique human genes (Supplementary Table 16). GO analysis revealed that these genes were enriched for biological processes that include regulation of the cell cycle, cell response to stress, and gene expression ( Fig. 2A). DGCR8-down-regulated miRNA gene targets were significantly enriched for genes expressed during fetal development across brain regions (Fig. 2B) and were not associated with any specific cell type (Fig. 2C), consistent with this biological pathway enrichment profile. Expanding analyses to to additionally include include gene targets of miRNAs down-regulated in hippocampus due to DGCR8 deficiency (Stark et al. 2008;Earls et al. 2012) yielded similar results (Fig. S4, Fig. S5).

Discussion
22q11DS is a rare genetic disorder characterized by neuroanatomic abnormalities that include widespread reductions in cortical SA and thicker cortex overall, with focal thinning in the caudal anterior cingulate, superior temporal cortex, and parahippocampus. As the typically deleted 22q11.2 region (LCR A-D) spans a gene-rich region, identifying which individual genes underlie these neuroanatomic alterations has remained a challenge. Here, by systematically examining the spatial convergence between severity of cortical SA and CT deviation in 22q11DS and expression of each 22q11.2 gene, we provide novel evidence prioritizing DGCR8 as a potential contributor to pervasive SA reductions in 22q11DS. AIFM3 was also associated with 22q11DS SA reductions. Interestingly, P2RX6 was associated with CT, suggesting that this gene may contribute to increases in CT in individuals with 22q11DS. Notably, regional expression of DGCR8 was robustly associated with severity of SA in 22q11DS across all age subgroups examined, suggesting that DGCR8 haploinsufficiency may disrupt early aspects of corticogenesis that are established by childhood.
DGCR8 is essential for the biogenesis of miRNAs, which regulate gene expression at the protein level (Rajman and Schratt 2017), and has previously been suggested to play a key role in 22q11DS phenotypes (Stark et al. 2008;Earls et al. 2012;Merico et al. 2014;Eom et al. 2020;Forsyth et al. 2020). Specifically, DGCR8 forms a complex with Drosha to cleave long primary miRNAs (pri-miRNAs) into short precursor miRNAs (∼70 nucleotides in length), before they are cleaved Figure 2. Characterization of 6804 unique gene targets of 59 miRNAs down-regulated in mouse cortex due to DGCR8 deficiency. A) Top five significantly enriched biological process, molecular function, and cellular component GO terms; B) enrichment for specific developmental periods and brain regions; and C) enrichment for specific CNS cell types, defined at varying specificity indices using the Specific Expression Analysis tool (Dougherty et al. 2010). Varying specificity thresholds in (B) and (C) are represented by the hexagon ring layers going from the least specific gene-lists (outer hexagons) to the most specific gene-lists (center), with hexagons scaled to the size of the gene-lists. BH-corrected Fisher's Exact p-values are plotted for each specificity threshold by color.
into mature miRNAs (∼22 nucleotides in length) by Dicer. Mature miRNAs are then loaded into the RNA-induced silencer complex (RISC), where they serve as a guide RNA strand for binding target mRNA transcripts to silence their translation into proteins. miRNA activity is increasingly recognized as a powerful posttranscriptional gene regulatory mechanism for diverse cellular processes, and in the brain, it is known to play a particularly important role in regulating cell proliferation, growth, and differentiation during fetal development (Yao et al. 2012;Rajman and Schratt 2017). Our finding that the gene targets of miRNAs down-regulated in mouse cortex due to DGCR8 deficiency are enriched for fetal-specific expression and regulation of the cell cycle and cell proliferation is consistent with these known roles of miRNAs. Notably, cell cycle parameters govern the balance between stem cell proliferation and differentiation during early brain development and thereby critically modulate cortical size and structure (Rakic 2009). In line with this, while knockout of DGCR8 is lethal, conditional knockout of DGCR8 in mouse embryonic stem (ES) cells has been found to disrupt cell cycle progression, ES cell proliferation, and differentiation (Wang et al. 2007). Similarly, knockout of DGCR8 in neural progenitors was found to disrupt progenitor pool maintenance, differentiation, and corticogenesis . Knockout of DGCR8 in pyramidal neurons also resulted in severe microcephaly, cell loss, altered inhibitory synaptic transmission, and premature death (Hsu et al. 2012). Although hemizygous depletion of DGCR8 yields milder effects on miRNA biogenesis relative to DGCR8 knockout, with some differences in brain phenotypes , DGCR8 +/− mice have been found to have altered cell proliferation, neurogenesis, and neuronal morphology (Stark et al. 2008;Fénelon et al. 2011;Ouchi et al. 2013;Amin et al. 2017), as well as enlarged brain ventricular volumes (Eom et al. 2020), altered homeostatic and synaptic plasticity (Earls et al. 2012;Amin et al. 2017), and cognitive deficits (Stark et al. 2008;Fénelon et al. 2011;Ouchi et al. 2013). Human induced pluripotent stem cell-derived neurons with hemizygous DGCR8 loss also show changes in calcium signaling and excitability (Khan et al. 2020). Our finding that DGCR8 expression is robustly associated with severity of SA deviation in 22q11DS across age subgroups is consistent with this literature implicating DGCR8 as a key contributor to brain and behavioral phenotypes in 22q11DS, and suggests that DGCR8 hemizygosity may contribute to SA reductions by disrupting miRNA modulation of cell cycle regulation during early brain development in 22q11DS.
AIFM3 was also associated with SA deviation in 22q11DS. AIFM3 is a proapoptotic protein that appears to stimulate cell death by depolarizing the membrane potential of mitochondria and activating the classical caspase-dependent apoptotic cascade (Xie et al. 2005;Zheng et al. 2019). Notably, in addition to the fundamental role of cell cycle regulation in defining the size and structure of the cortex, apoptosis also critically regulates cortex size (Haydar et al. 1999). Thus, two prominent waves of apoptosis are known to occur during corticogenesis. The first wave involves the elimination of a large number of dividing neuronal precursor cells during the peak of neurogenesis, regulating the size of the neuronal precursor pool. The second wave involves the elimination of postmitotic neurons during neuronal migration, regulating the wiring of developing neuronal networks (Blomgren et al. 2007). Although little is known about the specific role of AIFM3 in brain development, AIFM3 is highly expressed in the human brain (https://gtexportal.org/home/), supporting the possibility that AIFM3 haploinsufficiency could contribute to SA deficits by altering normal apoptotic processes during corticogenesis.
Finally, P2RX6 was associated with CT deviation in 22q11DS. P2RX6 encodes the P2X6 receptor-a member of the P2Xpurinergic family of ATP-gated ion channels that mediates fast excitatory postsynaptic potentials in neurons and smooth muscles (Motahari et al. 2019) and exerts neuromodulatory functions (Khakh and North 2012). While little is known about the role of P2RX6 during brain development, it is expressed in both developing and adult brain and is alternately spliced in the developing mouse brain and during in vitro neuronal differentiation (da Silva et al. 2007). Although speculative, P2RX6 deficiency could contribute to CT abnormalities in 22q11DS by altering modulation of synaptic signaling. Additional work is needed to experimentally validate that P2RX6 deficiency contributes to abnormalities in CT in 22q11DS and clarify underlying mechanisms.
Some potential limitations to the present study should be noted. First, our analyses are correlational in nature and assume causal, independent effects of deficient expression of one or more 22q11.2 genes on the regional severity patterns of SA and CT alterations in 22q11DS. Although the 22q11.2 deletion is known to yield significantly reduced expression of the majority of genes in the region (Stark et al. 2008;Jalbrzikowski et al. 2015;Lin et al. 2016;Gordon et al. 2019), experimental validation is necessary to confirm the present prioritization of DGCR8, AIFM3, and P2RX6 as genes that may be causally related to cortical alterations in 22q11DS. In addition, our analyses do not address the possibility that interactions between genes within the locus or with additional downstream interacting partners may also contribute to neuroanatomic alterations in 22q11DS, nor for the possibility that some brain abnormalities in 22q11DS may arise in part as downstream consequences of abnormalities in other organs, such as heart anomalies, which could be mediated by hemizygosity of genes within the 22q11.2 locus that are not expressed in brain (Schaer et al. 2009(Schaer et al. , 2010. In line with this, the strongest correlations (e.g., for spatial convergence between DGCR8 expression and SA) were moderate in strength, indicating that some variability in regional SA and CT deviance is not predicted by individual expression of these genes, as captured by the AHBA. Whether this is due to measurement error, interactions between genes and/or with downstream interacting partners, nonuniform effects of hemizygosity of individual 22q11.2 genes on corticogenesis in different brain regions, or other explanations remains to be determined. Additionally, we lacked high spatial resolution gene expression data from the developing human brain when SA and CT are initially shaped. We therefore relied on spatial expression data from the AHBA to derive insight into molecular mechanisms underlying neuroimaging phenotypes in 22q11DS, similar to other groups that have leveraged the AHBA to understand neuroanatomic abnormalities in neuropsychiatric, neurodevelopmental, and neurodegenerative disorders Romme et al. 2017;Fornito et al. 2019). Given these limitations, negative findings should be interpreted with caution, and future studies using gene expression data from developing brain samples collected at high spatial resolution will be important to confirm the present findings. Finally, our exploratory analysis of down-regulated miRNA gene targets due to DGCR8 hemizygosity was based on findings in postnatal cortical samples from 22q11.2-and DGCR8-deficient mouse models; it would be ideal to examine down-regulated miRNAs in cortex derived directly from individuals with 22q11DS or with DGCR8-specific hemizygous mutations, particularly during early corticogenesis. In the absence of such tissue, we nevertheless believe it is useful to identify biological pathways associated with gene targets of miRNAs previously identified as downregulated in cortex due to DGCR8 deficiency. Conversely, the large sample of individuals with confirmed A-D deletions and matched control subjects used to derive spatial measures of SA and CT deviation is a relative strength of the present study.
In summary, by integrating comprehensive maps of genomewide gene expression in the human brain and neuroanatomic data from the largest existing sample of 22q11DS individuals with molecularly confirmed A-D deletions, we prioritized DGCR8 and AIFM3 as potential contributors to cortical SA alterations in 22q11DS, and P2RX6 as a potential contributor to CT alterations. While DGCR8 deficiency has been found to modulate cell cycle progression, neural progenitor differentiation, and corticogenesis in animal models (Wang et al. 2007;Marinaro et al. 2017;Hoffmann et al. 2018), experimental validation is needed to confirm that AIFM3 and P2RX6 deficiency yield abnormalities in cortical SA and thickness, respectively. Nevertheless, DGCR8 and AIFM3 are involved in regulating two neurodevelopmental processes thought to fundamentally modulate brain size and structure (i.e., cell proliferation and apoptosis, respectively). Our finding that systematic investigation of 22q11.2 genes prioritized relatively understudied genes (i.e., AIFM3 and P2RX6) and a gene previously hypothesized to play a key role in 22q11.2 phenotypes (i.e., DGCR8), provides important opportunities to pursue novel mechanistic hypotheses in 22q11DS. Together, these results demonstrate the utility of combining neuroanatomic and publicly available transcriptomic datasets to derive mechanistic insights and prioritize individual genes that may underlie neuroanatomic differences in multigene CNVs.

Notes
Drs Ching and Thompson have received partial research support from Biogen, Inc. (Boston, USA) for work unrelated to the topic of this manuscript. Dr Antshel has received investigator-initiated research funds from Shire and has participated in an advisory panel for Arbor Pharmaceuticals. Dr D.G. Murphy has served on an advisory board for Roche. Dr Vila-Rodriguez has received research support from Brain Canada, the Canadian Institutes of Health Research, the Michael Smith Foundation for Health Research, and the Vancouver Coastal Health Research Institute and in-kind equipment support from MagVenture, and he has served on an advisory board for Janssen. Dr Linden receives editorial fees from Elsevier and book royalties from Springer Nature and Oxford University Press. Dr Owen has received research support from Takeda. Dr van den Bree has received research support from Takeda. Dr Thompson has received grant support from Biogen. The other authors report no financial relationships with commercial interests.