We investigated co-expression and regional patterns of gene expression in the healthy brain across the cortex, subcortex, and cerebellum using microarray gene expression data from the Allen Human Brain Atlas (AHBA)¹, a publicly available microarray dataset widely used for exploration of gene networks in the human brain (Hawrylycz et al., 2012). The microarray data was sampled from six adult human donors (3 White, 2 African American, 1 Hispanic, aged 24−57 years) in roughly 500 tissue samples from each donor, either in the left hemisphere only (n = 4) or in both hemispheres (n = 2). Although the same anatomical regions were sampled in all six donors, both the exact number and position of samples varied among donors. Donors also differed in other characteristics, including cause of death, post-mortem intervals, brain pH, tissue cytoarchitectural integrity, RNA quality, and number of probes used for each gene. Given space constraints, we refer the reader to the original technical white paper for additional details regarding the dissection methods, quality control, and normalization measures taken.²

The AHBA data were preprocessed and mapped to parcellated brain regions using a publicly available “abagen” processing pipeline.³ We applied the recommended default parameters outlined in the original manuscripts (Arnatkeviciute et al., 2019; Markello et al., 2021). Briefly, all available probes (Custom and Agilent) were included in the analyses. Probes that did not exceed background noise in at least 50% of all cortical and subcortical samples across all subjects were excluded. As more than one probe can be available for a single gene, the probe with the higher differential stability score was selected (default parameter). Probe selection was performed for each donor separately. The Montreal Neurological Institute (MNI) coordinates of tissue samples were updated to those generated via non-linear registration using Advanced Normalization Tools (ANTs). Tissue samples were assigned to brain regions in the provided atlas if their MNI coordinates were within 2 mm of a given parcel. To reduce potential misassignment, sample-to-region matching was constrained by hemisphere and gross anatomical divisions (cortex, subcortex/brainstem, and cerebellum). All tissue samples not assigned to a brain region in the provided atlas were discarded. We used the Desikan “aparcaseg” atlas (34 nodes per hemisphere + subcortex) to map tissue samples to cortical and subcortical regions⁴ (Desikan et al., 2006). We used the Diedrichsen atlas to map tissue samples to cerebellar regions⁵ (Diedrichsen et al., 2009). Inter-subject variation was addressed by normalizing tissue sample expression values across genes using a robust sigmoid function. Normalized expression values were then rescaled to the unit interval. Gene expression values across genes were normalized separately for each anatomic division (cortex, subcortex, and cerebellum) also using a robust sigmoid function. Samples assigned to the same region were averaged separately for each donor and then across donors. Gene expression values for the same region and gene sampled from both hemispheres from the same donor were averaged.

After implementing the pre-processing and quality control steps outlined above (Figure 1A), gene expression values from 15,633 protein-coding genes from the six donors were included in the analyses. Gene expression values were computed for 51 brain regions (34 cortical regions, 7 subcortical regions, and 10 cerebellar regions).

Median expression values (Log10 transcripts per million) for 10 distinct brain regions (amygdala, anterior cingulate cortex, caudate, cerebellar hemisphere, cerebellum, cortex, frontal cortex, hippocampus, nucleus accumbens, putamen, and substantia nigra) were obtained from the Genotype-Tissue Expression (GTEx) project database for the purpose of independent sample validation. Despite this dataset having lower spatial resolution compared to AHBA, we conducted a general comparison of co-expression patterns across the available regions for genes of interest to substantiate the generalizability of the AHBA expression data. We performed a Mantel test with permutations to examine whether the observed correlation between the two co-expression matrices was statistically significant. The permutation process involved shuffling the data multiple times to establish a null distribution, against which the actual correlation could be compared. P-values were calculated based on how often correlations in the permuted data are observed that are as extreme as, or more extreme than, the observed correlation.

Nineteen symptomatic C9orf72 HRE carriers participated in this study (10 males, 9 females; age range = 48−81 years, mean = 64 years, SD = 10 years). Individuals were recruited from the UCSF Memory and Aging Center (MAC) FTD cohort. C9orf72 HRE carriers were clinically diagnosed with FTD (n = 8), and FTD-ALS/ALS (n = 5), mild cognitive impairment (n = 4), and other (n = 2). For the two patients categorized as other, both showed concerns for motor neuron disease and possibly FTD. Healthy controls included 23 related family members (6 males, 17 females; age range = 28−72 years, mean = 47 years, SD = 13 years). Detailed information on participant inclusion criteria can be found in prior reports (Lee et al., 2014; Bonham et al., 2023). The UCSF Committee on Human Research approved the procedures for all participants. All participants or their surrogates provided informed consent prior to participation.

Genomic DNA was extracted from whole blood according to standard procedures. Participants were identified as carrying a pathogenic HRE in C9orf72 if they harbored >30 hexanucleotide repeats (DeJesus-Hernandez et al., 2011). Participants in this study were negative for pathogenic variants in MAPT and GRN.

All MRI scans were acquired on a 3T MRI scanner at the Neuroscience Imaging Center at UC San Francisco using previously described sequences (Bettcher et al., 2012). A high resolution T1-weighted image was acquired for all participants for structural reference, for the purpose of normalization, and for deriving morphological measures. MRI scans were processed using the FreeSurfer software package, version 6.0 (see text footnote 4) (Fischl, 2012). All images were visually inspected for segmentation accuracy and corrected as needed. We quantified disease burden using morphological measures of cortical thickness, since it provides a more sensitive measure of atrophy than gray matter volume (Winkler et al., 2010; Broce et al., 2023). For participants with multiple scans, the earliest scan with the best quality was selected. For visualization purposes, cortical brain images were generated using the R software statistical package “ggseg,” which allows for plotting brain atlases using simple features.

The R statistical package (version 4.2.2) was used for all statistical analysis. We first evaluated how average C9orf72 expression in each brain region varied across different anatomical divisions: cortex, subcortex, and the cerebellum. One sample t-tests (two-tailed) were conducted to assess which of the 51 brain regions from the six donor samples expressed mRNA to a significantly greater or lesser degree compared to average mRNA expression across the whole brain. To correct for multiple tests, reported p-values were Holm–Bonferroni adjusted. Cohen’s d values for one-sample t-tests were calculated to yield a measure of effect size.

We then evaluated C9orf72 co-expression patterns by correlating average C9orf72 expression values for each region with average expression values from 54 other genes known to cause ALS, FTD, or combined ALS-FTD (Table 1). These genes were selected based on prior reports (Abramzon et al., 2020; Kirola et al., 2022). Two ALS/FTD genes, PRPH (encoding Peripherin), and DAO (encoding D-amino-acid oxidase), were excluded from analyses at the preprocessing stage due to low quality or coverage. We separately correlated the columns of each expression matrix to generate a 54 × 54 gene co-expression matrix, reflecting the relative expression patterns across cortical, subcortical, and cerebellar regions for each gene pair (Figures 1B, C). We used the biweight midcorrelation as the similarity measure. In biweight midcorrelation gene expression, values that are much higher or much lower than the median value are given less weight in the correlation calculation, which makes the method more robust to extreme values and outliers (Langfelder and Horvath, 2008). Clusters were identified to assess co-expression patterns using the complete linkage method (Murtagh and Contreras, 2017).

Subsequently, to screen in an unbiased manner for novel genes that are co-expressed with C9orf72, we used the identical procedure described above. We correlated average C9orf72 expression values for each region with average expression of all available 15,633 protein-coding genes (Figures 1B–D). To explore whether C9orf72 co-expression patterns were driven by gene expression differences among the different anatomical divisions (cortex, subcortex, and cerebellum), we repeated the correlation analyses, separately, for each anatomical division. To reduce the influence of non-informative genes on subsequent analyses, genes with correlation values lower than an absolute value of 0.5 were excluded.

A main goal of this study was to identify a network of genes or gene products that may contribute to the regional vulnerability of the human brain to C9orf72 HRE-mediated pathology and, more generally, risk for ALS and FTD. Therefore, we identified correlations between average mRNA expression for each of the C9orf72-co-expressed genes from AHBA and differences in average cortical thickness from symptomatic C9orf72 HRE carriers compared to controls. Since age and sex are known to affect cortical thickness, and the age range within C9orf72 HRE carriers and controls was wide, our analyses controlled for age and sex. We followed established practices in the context of brain disorders to conduct these imaging transcriptomics analyses (Arnatkeviciute et al., 2022). In particular, our approach involved fitting separate linear models for each brain region to extract the estimated coefficients for case-control status on cortical thickness (i.e., C9orf72 HRE carriers versus controls), while also controlling for age and sex. We then correlated the coefficients obtained from the linear models for each region with gene expression across the same regions. To evaluate the significance of these correlations, we performed a permutation test with 10,000 iterations. This approach determines whether the observed imaging transcriptomics correlations differ significantly from what might occur due to chance.

We performed pathway enrichment analysis using the R statistical package “pathfindR” (Ulgen et al., 2019). Since enrichment analysis of a list of significant genes alone may not be informative enough to explain underlying disease mechanisms, we used “pathfindR,” which leverages interaction information from a protein-protein interaction network (PIN) to identify distinct active subnetworks and then perform enrichment analyses on these subnetworks. We preformed pathway enrichment analysis for three gene sets: “KEGG,” “GO-BP,” and “GO-MF” (all for Homo sapiens). For visualization, we plotted the top 30 terms for each gene set based on p-value.

Figure 2 displays how regional C9orf72 expression varies across the cortex, subcortex, and the cerebellum. The caudate was the only brain region that was significantly different from average expression across the whole brain after Holm–Bonferroni correction (p-adjusted = 0.006). The full results can be found in Supplementary Table 1. Group-level analysis tends to obscure how expression and imaging markers co-vary across regions within individual participants (Broce et al., 2023). In our subsequent analyses, we directly measure the strength of regional correlations at the individual-participant level, thus avoiding the dilution of these relationships that can occur at the group level. Specifically, we explored regional C9orf72 expression across cortical, subcortical, and cerebellar regions using correlation network analysis.

To explore whether the patterns of C9orf72 mRNA expression in brain resemble the patterns of mRNA expression of other well-established genes known to cause ALS, FTD, or combined ALS-FTD, we correlated the average C9orf72 mRNA expression values across the whole brain (51 brain regions) with average mRNA expression values across the whole brain for 54 ALS/FTD associated genes (Table 1). As shown in the dendrogram in Figure 3, we identified three main clusters: an orange cluster with 21 gene members, a green cluster with 16 gene members, and a blue cluster with 17 gene members. C9orf72 was in the orange cluster and grouped with 9 other gene members: ATXN2, TBK1, EWSR1, MFSD8, CTCF, TARDBP, TIA, ALS2, and GLE1. C9orf72 mRNA expression across the brain most strongly correlated with ATXN2 (bicor = 0.67) and TBK1 (bicor = 0.66) mRNA expression. Interestingly, as shown in the 54 × 54 gene co-expression matrix, the orange cluster was anti-correlated with several gene members from the green and blue clusters, including GRN, TREM2, TYROBP, OPTN, SOD1, ANG, CCNF, C21ORF2, and VCP, among others. Overall, the strongest positive correlations in the 54 × 54 gene correlation matrix were between TREM2 and TYROBP (bicor = 0.87), HNRNPA2B1 and FUS (bicor = 0.86), and TARDBP and TIA (bicor = 0.83). The strongest negative correlations were between GRN and TIA (bicor = −0.75), SOD1 and TIA (bicor = −0.74), and TARDBP, and GRN (bicor = −0.74). The full results can be found in Supplementary Table 2.

To independently validate the correlation patterns of the 54 specified genes associated with ALS/FTD, we compared gene co-expression in the AHBA dataset with the GTEx dataset using a Mantel test (see “Materials and methods”) (Supplementary Figure 1). The results indicated a statistically significant correlation between the gene co-expression patterns in the two datasets (bicor = 0.16, p = 0.0001) (Supplementary Figure 2). Notably, in the dendrograms generated from both datasets, C9orf72, TBK1, and ATXN2 formed a tight-knit cluster at the lowest of branch heights, suggesting greatest similarities of gene expression between these genes across both datasets. These results add a layer of confidence to the findings, reinforcing the conclusion that the observed similarity in co-expression is likely not a random occurrence.

By testing a large number of genes and correlations it is possible to create a complex dataset that is open to multiple interpretations. Therefore, to strengthen the connection between correlated gene expression and potential pathogenic roles, we explore whether genes associated with specific neurodegenerative diseases, such as SMCR8 and WDR41, known to form a complex with C9orf72 as part of normal cellular physiology, and TREM2, TYROBP, and GRN associated with leukoencephalopathy and FTD, show correlated gene expression. These analyses revealed intriguing patterns (Supplementary Figure 3): While C9orf72 negatively correlated with WDR41 (bicor = −0.52), this result did not achieve statistical significance following permutation analysis (p-value = 1.00). No significant correlation was found between C9orf72 and SMCR8. However, SMCR8 and WDR41 were positively correlated (bicor = 0.33, p-value = 0.005). Additionally, GRN, TREM2, and TYROBP displayed strong positive correlations among themselves (bicor ≥ 0.65, p-value = 0.0001). These relations were preserved in the larger AHBA 54 × 54 ALS/FTD correlation matrix (Figure 3), where GRN, TREM2, and TYROBP formed a closely-knit cluster, while SMCR8 and WDR41 closely clustered. Further, these relations were generally maintained in GTEx, with the exception of GRN, which integrated into a larger cluster alongside TARDBP, C9orf72, and other genes (Supplementary Figure 1).

We used the same data driven approach described above to screen for new genes that were co-expressed with C9orf72, beyond the 54 known ALS/FTD genes. Thus, we correlated the average C9orf72 mRNA expression values across the whole brain with average mRNA expression values for each of the 15,633 protein-coding genes. Then, to determine whether the whole brain co-expression patterns were driven by gene expression differences between the anatomical subdivisions, we visualized the similarities in whole brain co-expression and expression patterns among the different anatomical subdivisions (Figures 4A–D). Visualizing the correlation values between C9orf72 and a few gene members from the orange, green, and blue clusters (Figure 2) revealed that patterns across the whole brain most closely resemble the same relationships within cortical areas. Thus, for the imaging transcriptomic analyses described below, we focus on C9orf72-co-expressed genes, defined as genes with a correlation value higher than the absolute value of 0.5 across one or more anatomical subdivisions (n = 9,789 genes). The full results can be found in (Supplementary Table 3).

To identify a network of genes that may contribute to the regional vulnerability of the human brain to C9orf72 HRE-mediated pathology, we calculated correlations between average mRNA expression for each of the C9orf72-co-expressed genes from AHBA and differences in average cortical thickness from symptomatic C9orf72 HRE carriers compared to controls (Figures 1E–G). Both brain mRNA expression values from the AHBA dataset and brain imaging measures from the C9orf72 HRE carriers and controls were mapped to the FreeSurfer average cortical surface, allowing for these imaging transcriptomic correlations.

After permutation testing (p-value < 0.05), average mRNA expression from roughly 20% of all C9orf72-co-expressed genes (n = 2,120 genes) significantly correlated with average cortical thickness in symptomatic C9orf72 HRE carriers, including C9orf72 and 7 out of the 54 ALS/FTD genes: TARDBP, ATXN2, NEFH, SETX, PFN1, VAPB, and UNC13A (Figure 5). Other correlations within the top 10 most significant include PLEKHH3, RGS14, LRCH1, and PCDHB10, notable for their prior implication in dementia and aging (Supplementary Table 4).

We explored whether genes with significant C9orf72 case-control imaging transcriptomic correlations (n = 2,120 genes; henceforth, “C9orf72 imaging transcriptomic genes”) from the previous analyses were enriched for different cell populations in the brain. To do this, we used the publicly available Cell-type Specific Expression Analysis (CSEA) tools, which extracts data from the Brainspan database.⁶ These genes were enriched for several cell types or systems previously implicated in ALS and FTD (Figure 6) (BH adjusted p-value < 0.05), including layer 5b cells (Genc et al., 2017; Nana et al., 2019; Vatsavayai et al., 2019), cholinergic motor neurons in the spinal cord and brainstem (Casas et al., 2013), and dopamine type 1 and 2 receptor-positive (Drd1 + , Drd2 +) medium spiny neurons of striatum (Riku et al., 2016; Sobue et al., 2018).

We performed pathway enrichment analysis using the R statistical package “pathfindR” (Ulgen et al., 2019) to identify the biological pathways and mechanisms that the 2,120 C9orf72 imaging transcriptomic genes were enriched for (Figure 7, Supplementary Table 5). The most relevant KEGG biological pathways included autophagy; the Rap1, CAMP, and ErbB signaling pathways; and additional neurogenerative diseases (i.e., Parkinson’s disease and prion disease). Further, the most relevant GO biological pathway terms included mRNA splicing, endoplasmic reticulum to Golgi vesicle-mediated transport, protein ubiquitination, cellular response to DNA damage, regulation of miRNA transcription, and macroautophagy. Between KEGG pathways and GO biological pathway terms, over 25 genes were driving the autophagy enrichment, including C9orf72. Lastly, the most relevant GO molecular function terms included DNA binding transcription factor activity, RNA and microtubule binding, and protein ubiquitination. The full list of pathways and corresponding genes can be found in Supplementary Tables 5.