Comprehensive transcriptomic analysis indicates brain regional specific alterations in type 2 diabetes

Identification of DAGs in 13 brain regions

The major workflow of present study was shown in Figure 1. Control samples were matched with T2D cases to avoid bias resulting from cofounding factors. A total of 304 T2D samples and 608 matched control samples from thirteen distinct brain regions were then obtained. No significant differences were detected between matched samples for most of the covariates (Supplementary Figure 1 and Supplementary Table 1). Sample distribution after matching was plotted for each brain region, as shown in Supplementary Figure 2.

Figure 1. Workflow diagram of the study design. First, transcriptomic profiles of 13 human brain regions were derived from GTEx dataset. Second, differential analysis was performed to investigate regional specific changes. Distributions and functional annotations of DAGs were analyzed subsequently. Finally, Co-expression analysis was performed in caudate and hippocampus to study core modules and hub genes. The volcano plot and heatmap were generated using random sampling data of caudate transcriptome.

To identify DAGs, we performed generalized linear regression for genes against T2D status in each brain region using DESeq2 [20], with known and surrogate variables adjusted. At a 5% false discovery rate (FDR), less than 10 or even no DAGs were identified in most regions, indicating that T2D had no notable effect on transcriptome of the majority of human brain. Nevertheless, regulation of a considerable number of genes were disturbed in two regions, namely caudate (2939 DAGs) and hippocampus (486 DAGs) (Table 1). Among all the DAGs, only 178 genes were dysregulated in more than one tissue, and expression of all the ‘multi-hit’ DAGs responded to T2D status in the same direction except for 3 of them. This result indicated that different brain regions might associate with T2D in a region-specific manner. The details of DAGs in each brain region are provided in Supplementary Table 2. Control cases and T2D samples can be distinguished roughly in unsupervised clustering heatmap (P-value: 2.19eE-03 for caudate and 7.78E-07 for hippocampus, Fisher’s exact test) using identified DAGs (Figure 2A and 2B).

Figure 2. Distributions of DAGs. (A and B) Heat maps of DAGs (row) on samples (column) were shown for caudate (A) and hippocampus (B). (C) Gene biotype categories of DAGs according to GENCODE. X-axis shows the proportion to the total DAGs in each brain region. (D) Chromosomal distribution of DAGs. Y-axis shows the proportion to the total DAGs in each brain region. (E) The violin plot comparing the intra-chromosomal distances between DAGs in different brain regions. ^***P < 0.001; ^****P < 0.0001, Wilcoxon’s test. CCs, distance of DAGs within caudate group; HHs, distance of DAGs within hippocampus group; CHs, distance of DAGs between caudate and hippocampus group.

In view of the fact that caudate had the maximum sample size among all the 13 brain regions and the marked influence of sample size on differential analysis, we cannot evaluate the impact of T2D on diverse brain regions according to the number of DAGs directly. Therefore, we randomly selected samples with size ranging from 10 T2D samples versus 20 matched controls (the minimal size in 13 brain regions) to maximal number of each brain regions, with five T2D samples and ten matched controls added each time, and then bootstrapped 100 times. As expected, caudate got the largest number of DAGs at a 5% FDR level, followed by hippocampus, regardless of the sample size (Supplementary Figure 3). This result indicated that diverse brain regions were altered differentially in T2D. Furthermore, caudate and hippocampus might be the brain regions with the most significant transcriptomic changes in T2D.

Distributions and functional annotations of DAGs in various brain regions

To study the distributions of DAGs in various brain regions, we first categorized the DAGs according to gene biotypes in GENCODE [21]. Protein coding genes account for majority of DAGs (83.9% of caudate DAGs and 83.5% of hippocampus DAGs, Figure 2C). Whilst there was also a small subset of DAGs comprised of non-coding genes, including antisense RNA, pseudogene, lincRNA and others, implying the function of these non-coding genes in pathological process of diabetic brain.

Next, we investigated the chromosomal distribution patterns of DAGs, which showed that they were widespread across chromosomes (Figure 2D). More than 95% of DAGs located on autosomes and the overall distribution of caudate and hippocampus DAGs were similar. However, the proportion of DAGs located on small, gene-rich chromosomes (chr16-22) was slightly higher in hippocampus than in caudate.

As genes involved with the same disease process tend to locate adjacently [22], we then calculated the genomic distance of each pair of DAGs on the same chromosome (Figure 2E). Results showed that distances of DAGs within hippocampus group were the smallest compared with distances of DAGs within caudate group (P-value=8.95E-14) and between groups (P-value=1.11E-4). Surprisingly, distances of DAGs between groups were also smaller than distances of DAGs within caudate group (P-value = 1.82E-27). In general, above results suggested regional specific distribution of DAGs in various brain regions.

To obtain a functional overview of DAGs in different brain regions, we dissected functional annotations of the up- and down-regulated DAGs respectively. As a result, top Gene Ontology (GO) biological process terms enriched with down-regulated DAGs in caudate and hippocampus both related to synaptic functions (Supplementary Table 3), which have been reported to impaired in diabetic brain [23]. Nonetheless, functions of up-regulated DAGs tend to be more regional specific. In caudate, up-regulated DAGs enriched in terms like cilium movement and cilium organization, implying the hidden association of T2D and ciliopathy [24]. However, for up-regulated DAGs in hippocampus, terms related to muscle function and muscle tissue morphogenesis were enriched. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis revealed that DAGs in caudate also enriched in neurodegenerative diseases pathways such as Alzheimer disease, Huntington’s disease and Parkinson’s disease. In consideration of aforementioned association between T2D and neurodegenerative diseases [6–8], the caudate nucleus and related DAGs might play a vital role in these processes.

Co-expression analysis of caudate and hippocampus

Although numerous DAGs were identified in caudate and hippocampus, limited information could be provided by individual genes. Moreover, DAGs determined by differential expression analysis methods were biased to genes with large change in expression, while overlooking the coordination of gene expression. To integrate the expression difference into a systems-level context, we performed weighted gene co-expression network analysis (WGCNA) on caudate and hippocampus. We detected 14 modules in caudate (Supplementary Figure 4 and Supplementary Table 4) and 40 modules in hippocampus (Supplementary Figure 5 and Supplementary Table 5). The average size of caudate modules was larger than hippocampus (P-value: 9.43E-03, Wilcoxon test), which implied a stronger co-expression of transcriptomic profiles in caudate.

Among the 14 co-expression modules in caudate, eigengenes of 4 modules (cyan, green, blue and purple) were positively correlated (up-regulated) with T2D, while black and turquoise module were negatively correlated (down-regulated) with T2D (Figure 3A, P-value < 0.05, Student t test). No modules correlated with any of the potential confounding variables which have been regressed out. Gene significance (GS) varied dramatically across modules, which also indicated the different levels of correlation between co-expression modules and T2D (Figure 3B). To examine gene constitution of particular modules, we plotted GS against module membership (MM) for T2D-associated modules (Figure 3C). Significant correlations were detected in black, turquoise and blue modules, illustrating that genes significantly associated with T2D status are often also the important elements of these modules. This result confirmed the crucial role of these modules in type 2 diabetes process.

Figure 3. Co-expression modules in caudate. (A) Each row corresponds to a module eigengene and column indicate T2D status. The table were colored by correlation according to the legend. Each cell contains the corresponding correlation and P-value. (B) Each bar indicates the average of gene significance measure for all genes in a given module. (C) GS vs. MM plot for modules significantly correlated with T2D status. Each point corresponds to an individual gene within a given module, which was plotted by GS on the y-axis and MM on the x-axis. The regression line, correlation value and P-value were shown for each module.

In hippocampus, five of the 40 co-expression modules positively associated with T2D and another five modules negatively associated (Supplementary Figure 6A, P-value < 0.05, Student t test). The mean GS values across modules were shown in Supplementary Figure 6C. Among the 10 T2D-associated modules, significant correlations of GS and MM were detected in 7 modules (Supplementary Figure 6B).

Cell-type and regional specificity of co-expression modules

Disorder of regional specific gene networks usually results in impairment of regional specific functions. Hence, it is of great value to determine whether brain regional specific modules were disturbed in T2D, which might help unravel the underlying mechanisms of T2D-associated brain disorders. By testing the enrichment of distinct regional and cell-type specific expression markers [25], several modules in caudate were found to be significant at pSI 0.0001 (Figure 4). For instance, the green module, positively correlated with T2D, was significantly enriched in cortical astrocytes markers. The red module was significantly enriched in expression markers of oligodendrocyte of cortex and cerebellum, whilst eigengene of this module did not show prominent correlation with the T2D status. The most surprising result is the turquoise module, which was found to be striatum specific at every pSI (Supplementary Table 6). Cell-type markers of two major neuronal populations in striatum: Drd1-expressing and Drd2-expressing medium spiny neurons (D1- and D2-MSNs) were also enriched in turquoise module at every pSI. In consideration of the transcriptional association between turquoise module and T2D, this module could be determined as a striatum-specific module which was disturbed in T2D status.

Figure 4. Enrichment of brain regional and cell-type markers in caudate modules. (Top) Enrichment of caudate modules in markers of various neuronal and glial cell types. (Bottom) Same as above, but using markers for different brain regions. Asterisks indicate significant enrichment after Bonferroni adjustment. ^*P < 0.05; ^**P < 0.01; ^***P < 0.001; ^****P < 0.0001, Fisher’s Exact test. Cholin., cholinergic; seroton., serotonergic; cereb., cerebellum; UPBCs, unipolar brush cells; oligod., oligodendrocyte; MSNs, medium spiny neurons; hypocret., hypocretinergic.

However, in hippocampus, there were no T2D-associated modules enriched in any regional or cell-type markers at pSI 0.0001 (Supplementary Figure 7 and Supplementary Table 7). The darkolivegreen module was the only module enriched in hippocampus markers. Nevertheless, eigengene of this module had a weak correlation with T2D. Therefore, the darkolivegreen module may not account for the alteration of hippocampus in T2D. Next, we focused only on caudate to further study the correlation between T2D and these putative core modules.

Identification of modules genetically associated with T2D in caudate

Though genes in modules specific to caudate were dysregulated in T2D, it was unclear whether these impairments were only consequences of T2D or genetically associated with T2D. To identify the causal link of modules to T2D, we tested for their enrichment in T2D and height (as a negative control) associated single nucleotide polymorphisms (SNPs) from large genome-wide association study (GWAS) data sets.

In caudate, the turquoise module exhibited highly significant enrichment in T2D risk SNPs, while no such significant result was found in height-associated SNPs (Figure 5A, 5B and Supplementary Table 8). As the turquoise module was also highly enriched in striatum and MSNs markers, it was of great interest whether the T2D risk SNPs signal of this module was in connection with these specific markers. To explore this, we extract sub-modules of the turquoise module which only contain striatum, D1-MSNs or D2-MSNs specific markers at pSI 0.01 (sub-module derived at pSI 0.001 and 0.0001 was too small), respectively, and then performed the same enrichment method. All the three sub-modules were significantly enriched in T2D risk SNPs compared with height, especially for D2-MSNs. It can be inferred that the genetic risk for T2D associated with the turquoise module might distribute mainly across striatum and MSNs specific genes. Another 4 modules (red, magenta, green and tan module) were also relatively significant compared to negative control. However, only eigengene of green module was correlated with T2D and this module was enriched in cell-type markers of astrocytes of cortex.

Figure 5. Caudate modules enriched in T2D genetic signals. (A) For each module (or sub-module), the null distribution of T2D SNPs enrichment scores generated by 20,000 random permutations is shown. Real enrichment scores were depicted by red vertical lines. Modules were considered significant if FDR < 0.05. (B) For each module in (A), enrichment FDR for T2D SNPs are shown by histogram compared to height SNPs. Y-axis was log10 transformed and broken axis was used to show zero value. The red horizontal line marks the FDR threshold for significance, which is 0.05. FDR, false discovery rate.

Functional annotations and hub genes of caudate modules

Previous studies have shown that co-expressed genes tend to be functional related [26], we therefore studied the functional annotations of interested modules. Actually, although more than half of input genes were unassigned to specific co-expression modules, there were still large overlaps between DAGs and several modules. For instance, the blue module accounted for nearly 40% of up-regulated DAGs in caudate and a half of down-regulated DAGs were assigned to turquoise module (Supplementary Table 9). Nonetheless, the co-expression modules also harbored non-T2D-associated genes and had a greater power to delineate T2D-relevant transcriptional changes compared with DAGs. As an illustration, similar but more significant pathways were enriched for the blue module than up-regulated DAGs (Supplementary Table 10). The turquoise was also enriched in terms related to synaptic functions more significantly (Figure 6A), compared with down-regulated DAGs. In terms of other aforementioned T2D-associated modules, the cyan module associated with ribosome and the purple module is an immune module, revealing corresponding dysfunction in T2D brain. It was worth mentioning that the green module enriched in many metabolic processes and axonogenesis, in accordance with its enrichment in cortical astrocytes marker.

Figure 6. Functional annotations and hub genes of caudate modules. (A) (Left) Top 20 GO biological processes significantly enriched in turquoise module. (Right) Top 20 KEGG pathways significantly enriched in turquoise module. Numbers in the parenthesis indicate the numbers of genes associated with the respective terms. (B) Network plots showing top 500 connections in turquoise module; genes with most connections (hub genes) are shown in center. The size of each dot is proportional to log2 (number of connections for each gene).

Other than taking the turquoise module as a whole, we further investigated enriched pathways of striatum and MSNs sub-modules. All of these were enriched in pathways related to learning, memory and cognition (Supplementary Table 10). Taken together, the turquoise module was considered as a core functional module of caudate and its related synaptic impairments may contribute to the cognition decline in T2D.

Finally, we tried to identify highly connected hub genes for the turquoise module (Figure 6B). The top two hub genes were NSF (N-ethylmaleimide-sensitive factor, Entrez ID 4905) and ADD2 (adducin 2, Entrez ID 119), both of which were highly expressed in brain. NSF is essential for neurotransmitter release, together with other key factors of synaptic fusion machinery, such as SNAREs (soluble NSF attachment protein receptors) and SNAP (soluble NSF adaptor protein). They have been reported to play an important role in both diabetes and neurodegenerative disorders [27]. ADD2 encodes the beta subunit of adducins, which are a family of cytoskeletal proteins. Polymorphism of ADD1 gene (encoding the alpha subunit, Entrez ID 118) has been reported to associated with T2D [28], while little is known about the relationship between ADD2 and T2D. Thus, further studies are needed to explore the roles of these hub genes in T2D brain.

GTEx tissues and subjects

The GTEx project (v7, released in June 2017) provides expression data of 13 human brain regions from 752 post-mortem donors. Two of the brain regions were initially sampled from cerebellum and cortex preserved using the PAXgene tissue preservation system, and another 11 regions were subsequently sampled from frozen brains, including amygdala, anterior cingulate cortex, caudate, cerebellar hemisphere, frontal cortex, hippocampus, hypothalamus, nucleus accumbens, putamen, spinal cord and substantia nigra. Details regarding the sample collection, RNA sequencing and data processing are available at GTEx consortium paper [36].

To draw reliable conclusion, we only keep high-quality sequencing samples with RINs > 6.0. Cases with type 1 diabetes or unknown T2D status, and races other than black or white were excluded from this study.

Confounding factors including age [37], gender [38], race, BMI [39] and RIN [40] have been reported to be correlated with gene expression. To avoid bias, we used an optimal matching algorithm in R package MatchIt [41] to balance them between control and T2D groups, with optimal ratio of 2:1. The number of remaining matched samples of 13 brain regions ranged from 30 to 108.

Identification of regional specific DAGs

DESeq2 [20] (v1.22.2) was employed on all of the 13 brain regions to identify regional specific DAGs using raw read counts. Independent filtering of genes with low read counts was performed automatically by DESeq2 with alpha=0.05, and genes remained were referred to as ‘detectable genes’. To correct the known covariates as well as remove inferred hidden confounders in GTEx expression data, we employed “svaseq” function in sva R package [42] to identify 3 surrogate variables. The known confounding factors and surrogate variables were then added to the DESeq2 design, and negative binomial (NB) generalized linear regression model (GLM) was performed. Benjamini-Hochberg (BH) algorithm was used to adjust the Wald test P-values for multiple testing. Raw counts data were transformed to continuous, homoscedastic regularized logarithm transformed (r-log) values for further analysis.

Genomic distance and functional annotations

Genomic distance of each pair of DAGs on the same chromosome was calculated according to their genome coordinate positions. Wilcoxon rank-sum test was used to compare the distributions of intra- and inter-groups.

Functional enrichment analysis was performed using R package clusterProfiler [43] to identify significant GO biological process terms and KEGG terms. The FDR adjustment for P-value was made using Benjamin-Hochberg procedure and an FDR cutoff of 0.05 was used.

Weighted gene co-expression network analysis

To regress out uninterested sources of large variation, linear models containing age, gender, race, BMI, RIN and surrogate variables were fitted on the r-log data for all of the ‘detectable genes’ in each brain region, respectively. Then co-expression analysis was performed on residuals using WGCNA R package [44]. The soft-thresholding power were picked according to the scale-free topology criterion, and a singed gene network was constructed using blockwiseModules function in a single block with parameters mergeCutHeight = 0.15, minModSize = 40 and minKMEtostay = 0.7.

The module eigengene (ME) is defined as the first principal component of a given module, which can be considered as a representative of the gene expression profiles in a module. Module membership (MM) is a measure of gene-to-module membership by correlating its gene expression profile with the module eigengene of a given module. Gene significance (GS) was defined as the correlation between the gene and the T2D status.

Module graphs showing the top 500 connections were plotted using the iGraph [45] R package. Hub genes were defined as top 5% in number of connections.

Cell-type and regional specificity analysis

R package pSI [46], which provides lists of expression markers for diverse brain regions and cell types, was used to perform Fisher’s Exact test for regional or cell-type marker enrichment in co-expression modules at different specificity index thresholds (pSI 0.01, 0.001 and 0.0001; pSI 0.05 was deprecated for increased false positives). Regional markers were derived from Atlas of the developing Human Brain (www.brainspan.org), while cell-type markers were originally identified using translational profiling of genetically tagged cell lines purified from mouse brain [25].

GWAS set enrichment analysis

Enrichment for GWAS signal was conducted as previously described [47]. SNPs and their associated P-values were from large GWAS data sets of T2D [48] and height [49]. SNPs were assigned to genes if located within gene boundaries with additional 20kb on 5’ end and 10kb on 3’ end. The most significant SNP of each gene was selected, and then all of the genes were ranked with associated scores (-log₁₀ P-value) to calculate gene set enrichment score (ES) based on the Kolmogorov-Smirnov running-sum statistic using GSEA-P [50]. The null distribution of ES was generated by 20,000 random permutations of gene labels and associated scores. Enrichment scores were scaled by subtracting the mean and dividing by the standard deviation of permutation scores to correct for the gene set size, and empirical P-values were determined by the resulting z-scores.

Data and resource availability

The GTEx gene expression data and phenotype data were downloaded from dbGaP (http://www.ncbi.nlm.nih.gov/gap) under accession number phs000424.v7.p2. Statistical analysis was performed on R (v3.5.0). Scripts are available at https://github.com/ZedekiahZhou/T2D_Brain. All data generated or analyzed during this study are included in the published article (and its online supplementary files).