Gene Expression has Distinct Associations with Brain Structure and Function in Major Depressive Disorder

2.1. Case‐Control Brain Differences

Using the REST‐meta‐MDD dataset with 848 MDD patients (534 females and 314 males, age ranges: 18–65 years old) and 794 healthy controls (467 females and 327 males, age ranges: 18–64 years old), we compared the case‐control differences of three major functional brain measures (amplitude of low‐frequency fluctuation (ALFF), fractional ALFF (fALFF), and regional homogeneity (ReHo)) in 210 cortical regions. We found similar case‐control difference patterns across the brain for these three functional measures, generally with increased neuronal activity in lateral regions and decreased activity in medial regions, although different statistically significant brain functional abnormalities were observed after false discovery rate (FDR) multiple testing correction (Figure 1A, Figure S1, Supporting Information). MDD patients had significantly decreased ALFF in the left medial superior frontal cortex and significantly increased ALFF in inferior parietal cortex and temporal regions (Figure S1A, Supporting Information), whereas the MDD patients had significant decreased fALFF in the medial superior frontal and caudodorsal/pregenual cingulate cortex (Figure S1B, Supporting Information). In MDD patients, ReHo was significantly increased in the lateral prefrontal cortex, and decreased in the dorsal insula, precentral/postcentral/paracentral cortex, and caudodorsal cingulate cortex (Figure S1C, Supporting Information). The effect sizes for these functional differences ranged from Cohen's d −0.2 to 0.2, Moreover, brain‐wide case‐control differences for ALFF, fALFF, and ReHo were highly positively correlated with each other (correlation r _s ranged from 0.43 to 0.71) (Figure 1B, Figure S1D, Supporting Information). We therefore used principle component analysis (PCA) to extract the first principal component (PC) of effect sizes for functional measures, defined as funcPC1, which explained 71% variance of functional brain abnormalities in MDD (Figure 1C).

GMV showed smaller effect sizes (Cohen's d: −0.1–0.12), with decreased GMV in frontal and temporal cortex and increased GMV in occipital and parietal regions (Figure 1D), but none of these case‐control differences were significant.

To verify whether the observed case‐control differences were dependent on the brain parcellation scheme, we repeated the analyses based on the Desikan–Killiany (DK) atlas with 68 cortical regions. These results showed a comparable pattern of structural and functional differences (Figure S2, Supporting Information).

2.2. Relationships of Gene Expression with Case‐Control Differences

We used a whole brain transcriptomic dataset, AHBA gene expression matrix based on BN atlas (210 regions × 15 633 genes). PCA was used to extract the first PC from the matrix (genePC1), which explained 18.6% of the variance (Figure 2A). We then estimated the association of genePC1 with funcPC1 as well as Cohen's d for GMV and the functional measures. GenePC1 was negatively correlated with funcPC1 (r _s = −0.281, p = 3.64 × 10⁻⁵) (Figure 2B), and genePC1 consistently showed negative correlation with effect sizes for ALFF, fALFF, and ReHo (Figure S3, Supporting Information). In contrast, genePC1 was positively correlated with GMV differences (r _s = 0.367, p = 4.36 × 10⁻⁸) (Figure 2C). The opposite relationships of whole gene expression with funcPC1 and GMV differences indicates that the distribution of gene expression in the cortex is differentially related to functional and structural difference patterns in MDD. This pattern of opposite relationships of gene expression with structural and functional brain abnormities was replicated using data from the UK Biobank (Figures S4 and S5, Supporting Information).

Because females are nearly twice as likely to develop depression than males, we explored whether there were sex‐related differences in the observed results. We therefore split the dataset into a male group (314 MDD patients and 327 healthy controls) and a female group (534 MDD patients and 437 healthy controls), and investigated the relationships of genePC1 with brain abnormalities in both groups separately. Although no significant structural or functional case‐control differences were observed in the two groups (presumably due to the decrease in statistical power), the distribution of interregional differences in the male or female group were similar to those identified in whole dataset (Figure S6, Supporting Information). Moreover, in both groups, gene expression showed opposite relationships with structural and functional differences, which were consistent with the results in the whole dataset (Figures S7–S9, Supporting Information). These results indicate that the pattern of opposite relationships did not show differences between males and females.

Moreover, age ranged from 18 to 65 years old in REST‐meta‐MDD dataset. To test age influences on the observed relationships, we divided the dataset into an older group (≥32 years old, 420 MDD patients, and 433 healthy controls) and younger group (<32 years old, 428 MDD patients, and 361 healthy controls) (Figure S10, Supporting Information), based on a median split. In both groups, the interregional structural and functional case‐control differences were highly positively associated with those identified in the whole dataset (Figure S11, Supporting Information). In the older group, we found opposite relationships of genePC1 with funcPC1 and GMV differences, which were consistent with the results in the whole dataset (Figure S12, Supporting Information). In contrast, genePC1 showed weak positive correlations with funcPC1 and GMV differences in the younger group (Figure S13, Supporting Information). These findings indicate that the relationships may tend to emerge in a middle‐aged population.

We then identified the transcriptional correlates of brain abnormalities in MDD. Expression on 210 cortical regions of each gene were related to the funcPC1 and GMV differences, respectively (Tables S1 and S2, Supporting Information). After Bonferroni multiple testing correction, gene expression of 412 genes showed significant positive correlations and 415 genes showed significant negative correlations with funcPC1 (funcPC1+ genes and funcPC1−gene respectively), (Figure 2D, Figure S14A, Supporting Information). It indicates that higher gene expression in cortical regions is associated with increased and decreased neuronal activity in MDD for funcPC1+ and funcPC1− genes, respectively. For GMV, we identified 338 positively and 298 negatively related genes, which we termed GMV+ genes and GMV− genes, respectively, indicating higher gene expression in cortical regions with increased and decreased GMV in MDD, respectively (Figure 2E, Figure S14B, Supporting Information).

2.3. Relationships of Transcriptional Correlates with Gene Dysregulation

Following the approach of Gandal et al.,^{[ 20 ]} we obtained the meta‐analytic estimates of differential gene expression (DGE) in postmortem brain tissues of patients with major psychiatric disorders including MDD, schizophrenia (SCZ), autism spectrum disorder (ASD), and bipolar disorder (BP) (Table S3, Supporting Information). We first related in vivo transcriptional correlates of brain abnormalities to post‐mortem gene dysregulation in cortex of patients with MDD for genes. DGE values of genes for MDD were negatively correlated with their correlations with funcPC1, but not significantly correlated with transcriptional correlates of GMV differences (Figure S15, Supporting Information). However, DGE values of most genes were distributed around zero, which could cause bias for estimating correlations, so we applied bin‐based correlation analysis to increase the statistical power.^{[ 17 ]} When ≈10 000 genes were divided into 100 bins, the genes with similar transcriptional correlates of brain abnormalities were merged into a bin, which avoided the DGE values of too many genes gathering at a narrow range around zeros.^{[ 17 ]} We found a significant negative correlation with the transcriptional correlates of funcPC1, but a significant positive correlation with the transcriptional correlates of GMV differences (Figure 3A,E). The finding indicates that the transcriptional correlates of in vivo MDD‐related brain abnormalities from healthy AHBA donors captured the information of differential gene expression in ex vivo brain tissue from MDD patients.

Through bin‐based correlation analysis, we further observed the DGE values of ASD and BD were not significantly associated with the transcriptional correlates of funcPC1 (Figure 3B,C), whereas the DGE values of SCZ showed a significant positive correlation with the transcriptional correlates of funcPC1 (Figure 3D), which was opposite to the finding by MDD in Figure 3A. In addition, the DGE values of ASD, BD, and SCZ showed strong negative correlations with the transcriptional correlates of MDD‐related GMV differences (Figure 3F–H), which were also in opposite direction with the finding by MDD in Figure 3E. This indicated that the transcriptional correlates of MDD‐related brain abnormalities can capture the gene dysregulation of other psychiatric disorders, and the direction of relationships for MDD were opposite to ASD, BD, and SCZ. Comparable results were obtained at other bin numbers: 50, 150, and 200 gene bins (Table S4, Supporting Information).

2.4. Gene Set Enrichment

We first conducted the tissue‐specific expression analysis (TSEA).^{[ 21 ]} All of our identified genes were significantly enriched for brain tissue, especially funcPC1+ and GMV− genes (Figure 4A). We also conducted the cell‐type expression analysis based on the study by Seidlitz et al.^{[ 22 ]} The GMV− genes were specifically enriched for astrocytes, whereas the other three gene sets were significantly enriched for inhibitory and excitatory neurons (Figure 4B). Moreover, the funcPC1+ and funcPC1− genes also showed significant enrichment for microglia and endothelial cells, respectively.

In addition, we selected up‐ and downregulated genes for psychiatric disorders with the largest 500 positive and negative DGE values. We tested whether genes transcriptionally related to funcPC1 or GMV differences were enriched for genes that were highly upregulated and downregulated in psychiatric disorders (Figure S16, Supporting Information). The funcPC1+, funcPC1−, and GMV− genes were enriched for MDD downregulated genes, but not for MDD upregulated genes. In addition, widespread significant enrichments were identified for genes that were up‐ and downregulated in ASD, BP, and SCZ. Of them, GMV− genes showed the strongest enrichment for upregulated genes for BP and SCZ. Comparable results were observed using other thresholds (the largest 300–800 positive and negative DGE values) to define up‐ and downregulated genes for psychiatric disorders (Figure S17, Supporting Information).

We also conducted enrichment analysis for biological pathways including gene ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. The top 10 significant enrichment pathways are presented in Figure 5 . For funcPC1+ genes, we observed enriched pathways such as synaptic transmission (i.e., modulation of chemical synaptic transmission, synaptic signaling, and synapse organization) and transport processes (i.e., regulation of ion transport, and negative regulation of transport) (Figure 5A). For funcPC1− genes, significant enrichment was mainly found for signaling pathways: Calcium signaling pathway, steroid hormone mediated signaling pathway, and enzyme‐linked receptor protein signaling pathway (Figure 5B). In contrast, GMV− genes showed different enrichment terms, such as response to metal ion, neuroactive ligand‐receptor interaction, and retrograde endocannabinoid signaling (Figure 5C), whereas GMV+ genes were only enriched for three pathways: inorganic ion transmembrane transport, sensory organ morphogenesis, and skeletal system development.

Participants

In the main analysis, publicly available data with brain scans from the REST‐meta‐MDD project were used.^{[ 39 ]} The project totally included 1300 patients diagnosed with MDD and 1128 healthy controls. From them, 848 patients diagnosed with MDD and 794 healthy from 17 sites in China were extracted for group comparison analysis (Figure S18, Supporting Information). Age and sex did not have difference between the MDD and control group (age: t = −0.75, p = 0.94, sex: χ ² = 2.8, p = 0.094). Individuals with MDD were diagnosed as having current depressive status by experienced psychiatric physicians based on ICD10 or DSM‐IV. All participants provided written informed consent and the project was approved by local Institutional Review Boards (for details see Table S5, Supporting Information).

The REST‐meta‐MDD dataset was comprised of East Asian individuals from China. To assess whether the results can be generalized to other ethnic populations, the analyses were replicated in a European‐ancestry sample from the UK Biobank (802 depressed individuals and 802 controls) (Figure S19, Table S6, Supporting Information), matched for age, sex, intelligence, and educational attainment. More information about the sample selection, and phenotypic measures in the UK Biobank sample can be found in Supporting Information.

Imaging Procedures

Structural and functional MRI data were collected from 17 scanning sites using different 3T or 1.5T scanners, and different parameters were used for imaging acquisition (Table S5, Supporting Information). A MATLAB‐ and SPM‐based pipeline, DPARSF,^{[ 40 ]} was used to preprocess the resting‐state functional MRI (rsfMRI).^{[ 39a ]} Main procedures included removing the first 10 time points, slice‐timing correction, realignment, regressing out head motion effects using the Friston 24 parameter model,^{[ 41 ]} and normalization from individual native space to MNI space using DARTEL tool.^{[ 42 ]} After realignment, individual T1 weighted images were co‐registered to the mean functional image and then segmented into gray matter (GM) and white matter. GM images were also modulated and spatially normalized by the DARTEL tool.

After preprocessing, three main functional brain measures (ALFF, fALFF, and ReHo) were calculated using rsfMRI data. Next, all functional images were Z standardized and then smoothed with a 3 mm full‐width at half maximum (FWHM) Gaussian kernel. Similarly, GM images were also smoothed with the 3 mm FWHM Gaussian kernel. Finally, smoothed GM and three functional images were mapped to the BN atlas to extract GMV, ALFF, fALFF, and ReHo of 210 cortical regions, respectively. More details about the imaging acquisition and processing can be found in the study by Yan et al.^{[ 39a ]}

The UK Biobank data were analyzed using a series of imaging processing steps which can be found in Supporting Information.

Estimation of Regional Gene Expression

The AHBA microarray DNA data (http://human.brain‐map.org), comprising gene expression measures from a total of 3702 tissue samples by probes were used.^{[ 11 ]} The tissues were collected from six human donor brains 24 to 57 years of age with no known neuropsychiatric or neuropathological history (five males and one female). Of the donors, there were three Caucasian, two African‐American, and one Hispanic. Then, the abagen toolbox^{[ 43 ]} in Python was used to create a whole‐genome expression atlas by probe‐to‐gene and sample‐to‐region strategies using the following steps: 1) intensity‐based filtering of microarray probes to remove those that do not exceed background noise of 50%; 2) selection of a representative probe for each gene across each hemisphere; 3) matching of microarray samples to brain parcels from the BN atlas; and 4) normalization, and aggregation within parcels and across donors.^{[ 12 ]} This resulted in a 2D matrix that showed gene expression levels of 15 633 genes on 210 BN cortical regions. In addition, another cortical parcellation with 68 cortical regions, the DK atlas, was used to extract AHBA gene expression matrix for replication analysis.

Statistical Analysis

Two‐sample t‐tests were used to obtain the case‐control differences in GMV, ALFF, fALFF, and ReHo based on the BN atlas. T‐scores were then converted to Cohen's d effect sizes for ease of interpretation. The covariates, including age, sex, education, and head motion, were regressed out of the brain measures. The FDR correction was applied to correct for multiple comparisons using the Benjamini and Hochberg method with a threshold of p = 0.05. For a given brain measure, Cohen's d across 210 cortical regions was used to quantify structural or functional abnormalities in MDD. PCA was used to extract the first PC of effect sizes for the three functional measures, defined as funcPC1 that explained most variance of functional abnormalities. The first PC was also extracted from the gene expression matrix, defined as genePC1.

To test the relationship of whole gene expression of 15 633 genes for 210 cortical regions with the corresponding brain abnormalities, the correlation of genePC1 with GMV differences and funcPC1 was examined using the Spearman's rank correlation analysis.

To replicate the whole transcriptome‐neuroimaging relationships, identical analyses were conducted in the UK Biobank dataset (Supporting Information). In addition, possible effects of sex and age were tested on the transcriptome‐neuroimaging relationships (Supporting Information). Moreover, relationships between cortical expression distribution of individual genes and brain abnormalities were examined, indicating individual gene's transcriptional correlates of MDD‐related brain abnormalities. A Bonferroni multiple comparison correction was applied to select the genes significantly positively or negatively related to GMV differences, defined as GMV+ genes or GMV− genes. Similarly, genes related to funcPC1 were defined as funcPC1+ or funcPC1‐ genes.

Differential Gene Expression Analysis in Major Psychiatric Disorders

Following the approach of Gandal et al.,^{[ 20 ]} the DGE estimates in postmortem brain tissues of patients with major psychiatric disorders were obtained. In short, through linear mixed effects modeling, standardized beta coefficients were calculated as DGE values, which quantified the degree that a gene is up‐ or downregulated for a given psychiatric population. The tissues were sampled mainly from the prefrontal and anterior cingulate cortex of MDD patients, and the prefrontal cortex only for other psychiatric patients.

The ≈10 000 overlapping genes between AHBA and Gandal et al.^{[ 20 ]} datasets were obtained, and examined the associations between transcriptional correlates of brain abnormities in MDD and DGE values for psychiatric disorders for these common genes using Spearman's rank correlation analysis. It had been suggested that linking transcriptional correlates to DGE values directly for all genes could cause the bias for estimating correlation because of the noise,^{[ 17 ]} the bin‐based correlation analysis (Supporting Information) was further applied. In short, the common genes were ranked by their correlations with brain abnormalities (transcriptional correlates) in descending order. The ranked genes were then divided into 100 bins, which merged the genes with similar transcriptional correlates of brain abnormalities into a bin. For example, the genes with highest positive and negative correlations were divided into the first and 100th gene bin, respectively. The average transcriptional correlates and DGE values were calculated for each bin. For 100 gene bins, the relationships of average transcriptional correlates of in vivo brain abnormalities were examined with ex vivo DGE values for psychiatric disorders. Alternative numbers of bins (N = 50, 150, and 200) were also applied to test the stability across choices of bin numbers.

Enrichment Analysis

Gene set enrichment analysis was applied to identify functional annotations of the identified genes related to brain abnormalities in MDD (GMV+ genes, GMV− genes, funcPC1+ genes, and funcPC1− genes). In TSEA,^{[ 21 ]} a specificity index probability (p _SI = 0.05, permutation corrected) and Fisher's exact test with FDR‐BH correction were used to determine how likely our identified genes were to be specifically expressed in a given tissue. Next, cell‐type expression analysis were conducted based on the study by Seidlitz et al.^{[ 22 ]} They provided gene sets related to seven canonical cortical cell classes: inhibitory neurons, excitatory neurons, astrocytes, microglia, endothelial cells, oligodendrocyte precursors, and oligodendrocytes (Table S7, Supporting Information). The GSEApy package^{[ 44 ]} was used in Python to identify the cell types that genes related to brain abnormalities were enriched for. Moreover, the genes with the largest 300–800 positive and negative DGE values were selected as up‐ and downregulated gene sets for psychiatric disorders. Using the enrichment analysis, this work tested whether up‐ and downregulated genes for psychiatric disorders were expressed most in cortical regions that showed structural and functional abnormalities in MDD. In addition, the Enrichr tool available online (http://amp.pharm.mssm.edu/Enrichr) was used to identify enrichment pathways for the gene ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways with FDR‐BH correction.^{[ 44 ]}