Whole-transcriptome brain expression and exon-usage profiling in major depression and suicide: evidence for altered glial, endothelial and ATPase activity.

Brain gene expression profiling studies of suicide and depression using oligonucleotide microarrays have often failed to distinguish these two phenotypes. Moreover, next generation sequencing approaches are more accurate in quantifying gene expression and can detect alternative splicing. Using RNA-seq, we examined whole-exome gene and exon expression in non-psychiatric controls (CON, N=29), DSM-IV major depressive disorder suicides (MDD-S, N=21) and MDD non-suicides (MDD, N=9) in the dorsal lateral prefrontal cortex (Brodmann Area 9) of sudden death medication-free individuals post mortem. Using small RNA-seq, we also examined miRNA expression (nine samples per group). DeSeq2 identified 35 genes differentially expressed between groups and surviving adjustment for false discovery rate (adjusted P<0.1). In depression, altered genes include humanin-like-8 (MTRNRL8), interleukin-8 (IL8), and serpin peptidase inhibitor, clade H (SERPINH1) and chemokine ligand 4 (CCL4), while exploratory gene ontology (GO) analyses revealed lower expression of immune-related pathways such as chemokine receptor activity, chemotaxis and cytokine biosynthesis, and angiogenesis and vascular development in (adjusted P<0.1). Hypothesis-driven GO analysis suggests lower expression of genes involved in oligodendrocyte differentiation, regulation of glutamatergic neurotransmission, and oxytocin receptor expression in both suicide and depression, and provisional evidence for altered DNA-dependent ATPase expression in suicide only. DEXSEq analysis identified differential exon usage in ATPase, class II, type 9B (adjusted P<0.1) in depression. Differences in miRNA expression or structural gene variants were not detected. Results lend further support for models in which deficits in microglial, endothelial (blood-brain barrier), ATPase activity and astrocytic cell functions contribute to MDD and suicide, and identify putative pathways and mechanisms for further study in these disorders.

Introduction

Major depressive disorder (MDD) is common, with a 1-year prevalence of 10–15%. The lifetime mortality due to suicide in MDD reaches 15%. There is a strong association between suicide and psychopathologies such as MDD (1), but those at risk also have a diathesis with neurobiological underpinnings related to the predisposition to suicide (2). Traditional neurotransmitter system abnormalities explain only some of the neurobiology of suicide (3), and studies must separate the pathogenesis of MDD from the diathesis for suicidal behavior. Exploratory genome-wide brain gene expression profiling offers an unbiased way of seeking molecular phenotypes associated with MDD and suicide.

Genome-wide brain expression profiling of suicide and depression using oligonucleotide microarrays have identified potentially responsible genes and molecular pathways (3). Previous brain expression profiling studies of non-psychiatric sudden death controls (CON) compared with MDD found altered expression of synapse-related genes (4), genes implicated in the inhibition of cell proliferation (5), neuropeptides involved in stress response (6), and genes involved in oligodendroglial differentiation (7). Studies of suicide and MDD have found altered expression of GABAergic, glutamatergic and ATP biosynthesis related pathways (8–10). However, hybridization-based methods have limitations, including reliance upon existing knowledge about genome sequence, high background levels due to cross-hybridization (11,12), and a limited dynamic range of detection due to both higher background and signal saturation (13). Critically, hybridization-based oligonucleotide microarrays and protein assays are limited to detecting particular (protein coding) mRNA transcripts and protein outputs of a gene, yet it is now evident that most multi-exon genes undergo alternative splicing, which greatly increases the functional diversity of protein species (14). It has been shown that RNA editing of serotonin receptor 2C (5-HT2C) results in receptor isoforms that activate G protein less efficiently, and has been previously linked to depression and suicide (15). However, it is unknown whether more common forms of RNA processing such as alternative splicing play a role in these disorders.

Next generation sequence (NGS) technologies allow sequencing of all RNA transcripts produced by a given gene, and hence provide information about the whole exome (16). Unlike hybridization-based approaches, RNA-seq is not limited to detecting previously identified transcripts, but also novel transcripts including many tissue specific non-protein-coding RNAs. Another advantage is that RNA-seq has minimal background noise because the cDNA sequences can be unambiguously mapped to unique regions of the genome. RNA-seq has been shown to be highly accurate for quantifying expression levels, confirmed by using quantitative PCR (qPCR) (17) and spiked RNA controls of known concentration (18). Moreover, RNA-seq allows for the detection of splice junctions and can be used to measure differences in splicing events (i.e. differential exon usage, DEU) (19).

In the current study, we applied NGS whole-transcriptome profiling (RNA-seq) to identify genes, miRNA species, and molecular pathways that are altered (via differential expression or differential exon usage) in suicide and major depressive disorder (MDD) in dorsolateral prefrontal cortex (BA 9). The dorsolateral prefrontal cortex is involved in regulation of impulsivity, deicision-making, cognitive control of mood and other executive functions related to suicidal behavior (20,21), and functional imaging studies have consistently found deficits in perfusion and cortical thickness in prefrontal cortex of suicides prior to death and suicide attempters (22,23). Whereas previous studies compared non-MDD suicides, sudden death non-psychiatric controls (CON), and MDD suicides (MDD-S) (24), here we compare MDD-S to MDD non-suicides (MDD) and CON. This allowed separate evaluation of effects of MDD (vs. CON) on gene expression by comparing MDD combined MDD-S vs. CON, and the effect of suicide by comparing MDD-S vs. CON combined with MDD. DSM-IV diagnosis was done by psychological autopsy using a validated, structured clinical interview for axis I and axis II disorders and all subjects were free of neuropathology and had negative toxicological screens for psychotropic medication and illicit drugs. In addition to exploratory, whole-genome pathway analysis (FDR corrected at p<0.1), we also applied a hypothesis-driven approach at a more lenient threshold (p<0.05 uncorrected) to examine specific pathways previously reported by a related study of MDD-S, CON and non-depressed suicides using oligonucleotide microarrays (8).

Materials and Methods

Samples

Fifty-nine clinical samples were obtained from the brain collection of The Division of Molecular Imaging and Neuropathology, at the New York State Psychiatric Institute and Columbia University. All procedures for brain collection and psychological autopsy were approved by the applicable Institutional Review Boards. Psychiatric diagnosis in the suicides and depressed individuals as well as absence of diagnoses in the controls was determined by the Structured Clinical Interview for DSM IV (SCID-I and II) as part of a validated psychological autopsy method described elsewhere (25). All subjects were selected because they died suddenly to avoid metabolic complications related to agonal effects. All brains were free of gross neuropathology and had negative brain toxicology for psychotropic, illicit psychoactive drugs and neurotoxic drugs. There were no diagnoses of Alcohol or Drug Use Disorders. Antemortem medication history for three months ruled out recent exposure to psychotropic medication and confirmed results of comprehensive peripheral and brain toxicology. Brain samples were dissected from Brodmann Area 9 as previously reported (26) with an attempt to enrich for gray matter in when dissecting tissue from BA9. However, samples inevitably also contained a small amount of white matter.

We applied whole-exome, gene-level analysis of count data using DESeq2 (27) to examine differential expression between 21 MDD-S (subjects with major depressive disorder and suicide), 9 MDD (subjects with MDD and no suicide), and 29 sudden death healthy control (CON) subjects with no MDD and no suicide (59 samples total). Small RNA analysis (miRNA differential expression) used an age and sex matched subset that included 9 MDD-S, 9 MDD and 9 CON. The main results reported throughout the text derive from the full dataset of 59 samples, and adjust for age and gender. Unfortunately RIN scores from two samples (one MDD-S and one CON) were not estimated due to an unexpected ribosomal ratio, low 18S or low RNA concentration. Given that RIN score captures the integrity of ribosomal RNAs but fails to measure total mRNA integrity, and is not necessarily a determining factor in generating good quality RNA sequencing data (personal communication with Peter Nagy, Director of Columbia University Genomics Core in the Department of Pathology), in the main text we present results in the full sample (N=59) without covarying for RIN score, and compared these results to an additional analysis that includes RIN score as a covariate in a subset of 57 samples (results listed in Supplement). See Supplementary Material for further discussion on accounting for the possible effects of pH, PMI and RIN. The raw RNA-seq data will be deposited in publically available and structured repositories (i.e. Gene Expression Omnibus, see Supplementary Materials for more details).

RNA sample preparation and sequencing

RNA was extracted from dorsal prefrontal cortex (BA9) as previously described (26,28) and according to the guidelines recommended by the NIH Roadmap Epigenomics Mapping Consortium (REMC). We used the Ambion’s mirVana miRNA isolation kit (#. AM1560, Life Technologies, Carlsbad, CA, USA) to isolate total RNA from postmortem brain tissue (Brodmann Area 9). The coding RNA library was generated using the TruSeq Stranded Total RNA Sample Prep kit (Illumina, San Diego, CA, USA) which includes rRNA depletion and chemical fragmentation. The miRNA libraries were generated using the TruSeq Small RNA Sample Prep Kit (Illumina, San Diego, CA, USA) following the manufacturer’s recommendations (no rRNA depletion or chemical fragmentation was performed).

Paired-end, strand specific sequencing for total RNA was performed on Illumina HiSeq 2500 with 100 bp read lengths, while single-end sequencing was performed for microRNA on Illumina MiSeq with 50 bp read lengths. For each clinical sample, raw RNA-seq reads were aligned and mapped to the Ensembl GRCh37 human reference genome and assembled using Tophat v2.0.9 resulting in BAM files for each of 59 samples. Between 17,000,000 and 57,000,000 reads were obtained for each sample, and ~75–90% of reads were successfully mapped to the genome for each sample. Of these reads, ~10–15% aligned to multiple genomic loci. Read statistics for each sample (subject) are listed in Supplementary Table 1. Note that at the time of writing of this manuscript, the most recent GTF file compatible with Tophat was GRCh37. Hence some of the IDs listed in the main text tables have been deprecated and are therefore not included (i.e. they been replaced by one or more IDs) in the most recent (GRCh38) EnsEMBL database. The archived GRCh37 database can be accessed at http://grch37.ensembl.org/index.html.

For microRNA, the resulting FastQ files (adaptor sequences already trimmed) were first converted to FASTA format and processed according to the Small RNA Analysis pipeline in CLC Genomics Workbench version 7.5 (http://www.clcbio.com/products/clc-genomics-workbench/ ) using default settings (i.e. mature length variants were allowed 2 additional or missing upstream and downstream base-pairs and 2 mismatched base-pairs, and alignment was strand-specific). Unique reads were extracted and counted and subsequently merged and annotated according to known miRNA species included in miRBase Release 21. Total counts for each miRNA (1,353 unique annotations with positive counts for one or more samples) were divided by the total library size to arrive at an expression value for each sample. As the miRNA data are count data which are not normally distributed, non-parametric statistical tests were applied. These expression values were submitted to one-way Kruskal-Wallis analysis of variance, a non-parametric version of 1-way ANOVA, with three levels: MDD-S, MDD and CON. The samples for miRNA group comparisons were matched for age and sex (Table 1), and RIN scores and brain pH also did not differ between the MDD-S, MDD and CON groups (N=9 each group, p-vals > 0.75, Table 1). Therefore these variables were not included as nuisance covariates in the above model.

Table 1

Demographics of Study Groups

Total RNA	MDD-S	CON	MDD
Sample size	N = 21	N = 29	N = 9	Statistic
Age (yrs, mean, sd)	52.0 (21.7)	43.5 (21.3)	57.7 (15.5)	F(2,56)= 2.0, p = 0.14
Sex (Number, % Female)	8 (38%)	6 (21%)	3 (33%)	ChiSq = 1.9, p = 0.39
RIN (mean, sd)	6.6 (1.5)	7.1 (1.3)	7.4 (1.1)	F(2,53)=1.2, p=0.32
PMI (hrs, mean, sd)	16.1 (7.0)	13.2 (4.7)	15.1 (4.4)	F(2,54)=1.6, p=0.22
Brain pH (mean, sd)	6.4 (0.3)	6.5 (0.3)	6.4 (0.3)	F(2,55)=0.8, p=0.46
Micro RNA
Sample size	N = 9	N = 9	N = 9
Age (Mean Years, sd)	59.2 (14.2)	56.7 (16.0)	57.7 (15.5)	F(2,24)= 0.06, p = 0.94
Gender (Number, % Female)	3 (33%)	3 (33%)	3 (33%)	ChiSq = 0, p = 1
RIN (mean, sd)	7.7 (0.9)	7.4 (1.3)	7.4 (1.1)	F(2,23)=0.3, p=0.77
PMI (mean, sd)	13.3 (6.0)	16.7 (3.8)	15.1 (4.4)	F(2,24)=1.1, p=0.36
Brain pH (mean, sd)	6.3 (0.3)	6.4 (0.4)	6.4 (0.3)	F(2,24)=0.2, p=0.82

MDD-S = depressed suicides, CON= sudden death non-psychiatric controls, MDD = depressed non-suicides; RIN=RNA integrity score; PMI=postmortem interval. Last column shows statistic when testing for group differences in the corresponding demographic variable in each row.

Gene-level differential expression analyses

We used DESeq2 v.1.6.3 (27) to assess differential expression between groups. DESeq2 ignores ambiguously mapped reads and applies a Generalized Linear Model (GLM) to count data assuming a negative binomial distribution. This approach allows for adjustment of nuisance effects (here age, gender and RIN score) and uses shrinkage estimation for dispersion (i.e. variance within in group) and fold changes (pooling information both across samples and across all genes) to improve stability and interpretation of estimates. See Supplementary Materials for R session information and primary R code used to generate these results.

For each BAM file resulting from Tophat2, read counts per gene were summarized using the ‘summarizeOverlaps’ function in the GenomicAlignments R library (29) and a transcript database derived from the GRCh37 human genome assembly. For each gene, read counts were fit to a GLM with log link function (~ age + sex + condition or ~age + sex + RIN + condition), where condition was the factor of interest with three levels: MDD-S, MDD and CON. A likelihood ratio test (LRT) was performed, whereby the likelihood of the above model given the data was compared to a reduced model which excluded the “condition” factor. The LRT was performed in order to test all three levels at once (equivalent to an omnibus F-test). For significant testing in the LRT, an analysis-of-deviance between the full model (with the condition factor) and the reduced model is computed by means of a χ2 likelihood-ratio test (19).

In addition to the above model in which each group was modeled separately, two additional models were estimated whereby by the three levels (MDD-S, MDD and CON) were reduced to two levels: one in which MDD-S and MDD groups were combined in order to estimate effects of depression (MDD-S + MDD vs. CON), and another in which MDD and CON were combined to estimate effects of suicide (MDD-S vs. MDD+CON). Here log2 fold changes (LFC, or log2FC) of between-group differences were estimated via contrasts of individual coefficients and associated p-values were generated. Note that pairwise group (i.e. MDD-S vs. CON, MDD-S vs. MDD, MDD vs. CON) and suicide (MDD-S vs. CON+MDD) and depression (MDD-S + MDD vs. CON) comparisons were contrasts of parameter estimates generated from the above models which included age and gender as nuisance covariates. Hence these contrasts are controlling for age and sex effects.

In DESeq2, significance testing for group comparisons uses the Wald test, whereby shrunken estimates of LFC are divided by its standard error, resulting in a z-score, which is then converted to a p-value. See (27) for more details. For both LRT and contrasts, p-values for genes surviving an independent filtering step based on overall mean read count were adjusted for multiple comparisons correction using Benjamini-Hochberg correction for false discovery rate (FDR) at threshold of p<0.1 (as per all defaults in DESeq2, see (27)). For main exploratory analyses we did not require a pre-specified effect size (i.e. a minimum LFC threshold was not applied). However LFC is reported for each gene surviving multiple comparisons correction pairwise group comparisons and pairwise group comparisons applied a minimum LFC threshold (10% fold change).

qPCR Validation

RNA Extraction

RNA was extracted from samples using the Qiagen RNeasy Mini Kit (Qiagen, Germantown, MD; Cat. no. 74104) according to the manufacturer’s protocol and included optional on-column DNase treatment. cDNA Synthesis: RNA was quantitated by spectrophotometer and subsequently, 1μg total RNA was used as template to synthesize cDNA with the High Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Foster City, CA; Cat. no. 4368814) with a reaction volume of 50μL. Gene Expression Assays: Gene expression assays were purchased from Applied Biosystems (Cat. no. 4369016) and are supplied as a 20x concentration mix (final primer concentration 900nM each; probe concentration 250nM). Specific assay numbers are as follows: ACTB: Hs01060665_g1, MTRNR2L8: Hs04235260_sH, CCL4: Hs04421399_gH, MRPS6: Hs01122386_m1, IL8: Hs00174103_m1, SERPINH1: Hs01060397_g1. Real Time PCR Analysis: Real time PCR was performed on the ABI 7900HT Fast Real Time PCR System (Applied Biosystems) using assays specific for each gene of interest. Each reaction well contained μL of TaqMan® Universal Master Mix II (Applied Biosystems, Cat#4440039), cDNA equivalent to 20ng of total RNA and 0.5μL of each specific gene expression assay in a reaction volume of 10μL. Cycling conditions were as follows: 95°C for 10 minutes for polymerase activation, followed by 40 cycles of 95°C for 15 seconds and 60°C for 1 minute. Data analysis was performed using Sequence Detection System software from Applied Biosystems, version 2.4. The experimental Ct (cycle threshold) was calibrated against the endogenous control product ACTB. Samples were analyzed for relative gene expression by the ΔΔCt method (30).

Gene ontology enrichment analyses

P-values from the LRT test for each gene were used in analyses of gene ontology (GO) using the gene score resampling (GSR) method in ErmineJ (v3.0.2) (http://chibi.ubc.ca/ermineJ) (31) with default values and 200,000 iterations for P-value estimation (main text). Gene Score Resampling (GSR) uses the full list of genes and their scores (i.e the unadjusted p-values obtained from the DESeq2 analysis). Thus an advantage of this approach is that genes that do not meet a statistical threshold for selection can contribute to the score. In addition, more information contained in the gene scores is preserved than in a standard overrepresentation analysis (ORA), because ORA is essentially rank-based, whereas GSR uses the gene scores themselves. See http://erminej.chibi.ubc.ca/help/tutorials/running-an-analysis-resampling/ for more details. Left and right tailed p-values were also computed from the “two-group” models and entered into ErmineJ in order to assess both up- and down-regulated GO terms in comparisons testing for effects of suicide (MDD-S vs. MDD+CON) and depression (MDD+MDD-S vs. CON). All P-values resulting from ErmineJ are two-tailed and are corrected for gene multifunctionality (32) and multiple comparisons (FDR < 0.1). Hypothesis-driven gene ontology analyses which examine specific pathways of interest based on previous studies (8) report uncorrected p-values (Supplementary Figure 1).

Additional correlation (CORR) analyses were conducted using ErmineJ in order to identify potentially spurious results explained by intercorrelated gene sets (as opposed to shared effects of group). For these analyses, the rlog transform in DESeq2 was used to produce expression values for each gene in each sample (27). The point of the rlog is to remove the dependence of the variance on the mean, particularly the high variance of the logarithm of count data when the mean is low. Rlog uses the experiment-wide trend of variance over mean, in order to transform the data to remove the experiment-wide trend. The resulting expression data matrix was input to ErmineJ for CORR analyses. The uncorrected p-value testing for correlation in genes within the same GO category are listed in the “Correlated Pval” column in main result table (main text only). The GO terms with a low p-value, positive findings (i.e. significant effect of group status) for these pathways could have resulted from spurious correlations, thus they are noted as such in the tables but not discussed as part of the main findings in the text.

Differential exon usage analysis

We applied DEXSeq v1.12.2 (19), an “exon-centric” analysis which explicitly tests for differential exon usage by comparing a fully specified GLM (which models an exonXcondition interaction) to a reduced GLM (which excludes the exonXcondition term) using a likelihood ratio test (LRT, see methods for DESeq2 analysis). For each gene, the read counts that map to unique splicing events (i.e. exon deletions, insertions or alternate start/end sites, here referred to as “exon bin”) are modeled as ~ sample + exon + exon:condition + age:exon + sex:exon in the full model, and ~ sample + exon + age:exon + sex:exon in the reduced model. Note that the “sample” term models the overall gene expression value for each sample (i.e. read counts mapping to the whole gene), while “exon” models the deviation about this mean for each exon bin (i.e. unique splicing event) in the gene. Age:exon and sex:exon co-vary for interaction effects of age and gender with “exon”. Age and sex were not included as main effects since overall gene expression is modeled separately for each individual sample, thus absorbing any variance attributable to inter-individual differences in age and sex. Thus this approach explicitly tests for evidence of differential exon usage between groups (here modeled using the exon:condition interaction term) after accounting for group differences in overall gene expression.

RNA-seq variant discovery

The best practices workflow for calling variants from RNA-seq data using the Genome Analysis Toolkit (GATK) v3.4, considered a gold standard for variant calling from NGS data (33), was applied to the Tophat realigned BAM files from each sample. Briefly, following Picard preprocessing steps (adding read group information, sorting, marking duplicates and indexing), the GATK tool SplitNCigarReads was used to split reads into exon segments (getting rid of Ns but maintaining grouping information), hard-clip any sequences overhanging into the intronic regions, and reassign mapping qualities. This step was followed by indel realignment and base recalibration (BSQR). Finally, insertion-deletions and SNPs for each sample were called using GATK (34,35) HaplotypeCaller tool, merged and then filtered based on Fisher Strand values (FS > 30.0) and Qual By Depth values (QD < 2.0). The final merged vcf file was annotated for predicted effects using SnpEff v4.1k (36) to classify variants such as missense, nonsense, splice site, synonymous, intronic, or stop gain/loss etc. and effected genes. SnpSift v4.1 (37) was then used to annotate variant calls with dbSNP ids (if available) and estimate case-control (i.e. SuivsNonSui or MDDvsNonMDD) differences in genotype frequencies for various inheritance models (i.e. dominant, allelic) using Fishers Exact test or Cochran-Armitage test. Variants for the top 34 genes were thresholded at p<0.05 uncorrected (with either dominant or allelic model), while whole-exome results were thresholded at p<0.0001 uncorrected. Multiple comparisons correction was also applied using False Discovery Rate (FDR) correction. For more details, see the available code in Supplementary Material and (http://gatkforums.broadinstitute.org/discussion/3891/calling-variants-in-rnaseq and http://snpeff.sourceforge.net/SnpSift.html). For a complete list of predicted effects and variant types according to Sequence Ontology see http://snpeff.sourceforge.net/VCFannotationformat_v1.0.pdf.

Results

Demographics

29 sudden death non-psychiatric controls (CON), 21 MDD-S and 9 MDD samples were included in the analyses. Mean age, sex ratios and RNA quality scores were not different between groups (Table 1). Effects of age and sex on transcript expression were included as nuisance covariates and adjusted for prior to subsequent analyses (see below). Although mean RIN score was above 6.5 in each group, additional analyses included RIN score (an RNA quality measure) as a covariate in order to ensure results were not confounded by differences in RNA quality.

Whole-transcriptome brain expression differences in suicide and depression

An ‘omnibus’ test for group differences (likelihood ratio test, LRT, see methods) identified 34 genes that survived correction for multiple comparisons (p<0.1 corrected, Table 2). Plots for the top four differentially expressed genes are shown in Figure 1. These include Humanin-like 8 (MTRNR2L8, also known as HN8) which was higher in depression (MDD+MDD-S vs. CON, log2FC = 0.67, adjusted p=8.2E-5) and higher in suicide (MDD-S vs. CON+MDD, log2FC = 0.72, adjusted p=4.1E-9); serpin peptidase inhibitor, clade H (heat shock protein 47), member 1 (SERPINH1) which was lower in depression (log2FC=-0.72, adjusted p=2.5E-5) and lower in suicide (log2FC=-0.40); interleukin 8 (IL8, also known as CXCL8) which was lower in depression (log2FC=-0.53, adjusted p=0.0006), and chemokine (C-C motif) ligand 4 (CCL4) which was lower in depression (log2FC=-0.59, adjusted p=0.0003). Two of these top four genes showed group differences in an analysis using age- and sex-matched subsample (9 MDD-S, 9 MDD and 9 CON, MTRNR2L8: LRT χ2 = 10.9, p=0.004, IL8: χ2 = 8.23, p=0.015, CCL4: χ2 = 1.57, p=0.45, SERPINH1: χ2 = 3.18, p=0.20, unadjusted p-values, data not shown). Exploratory results showing effect sizes (fold change) for pairwise contrasts between each group are listed in Supplementary Table 2 (minimum 10% fold change, p<0.001 uncorrected). In the likelihood ratio test (LRT) analysis that covaried for RIN score, five of the top 10 differentially expressed genes no longer survived p<0.1 adjusted, while 9/10 survived p<0.15 adjusted and all 34 genes survived p<=0.001 uncorrected (Supplementary Table 3). Genes that did not survive p<0.1 adjusted when covarying for RIN score are denoted in Table 2 with an ǂ symbol.

Table 2

Genes differentially expressed among groups (CON, MDD and MDD-S, Likelihood ratio test (LRT) for ~ age + sex + condition vs. ~ age + sex, p<0.1 FDR corrected). pvalue=uncorrected pvalue; padj= FDR corrected pvalue; In rightmost columns, SUI>nonSui=MDD-S > (MDD+CON); MDD>nonMDD=(MDD+MDD-S)>CON.

Ensembl ID	Gene Symbol	Mean Normalized counts			LRT (χ2)	pvalue	padj	SUI > NonSui		MDD > NonMDD
		CON	MDD	MDD-S				Wald stat	pvalue	Wald stat	pvalue
ENSG00000255823	MTRNR2L8	17.5	30.4	46.5	43.1	4.3E-10	1.5E-05	7.4	1.6E-13	5.9	3.2E-09
ENSG00000149257	SERPINH1	154.5	104.3	72.3	30.2	2.8E-07	1.0E-02	-4.1	4.1E-05	-6.0	1.8E-09
ENSG00000225217	HSPA7*	106.0	80.7	57.2	28.1	7.8E-07	1.4E-02	-3.6	3.7E-04	-5.5	3.1E-08
ENSG00000169429	IL8ǂ	69.8	49.7	36.6	24.9	3.9E-06	4.6E-02	-3.4	7.2E-04	-5.1	3.1E-07
ENSG00000129277	CCL4ǂ	15.4	10.9	7.8	23.0	1.0E-05	6.0E-02	-3.7	2.3E-04	-5.4	6.3E-08
ENSG00000243927	MRPS6	181.0	152.9	143.8	23.1	9.6E-06	6.0E-02	-1.8	6.9E-02	-4.4	1.1E-05
ENSG00000259079	RP1-261D10.1***	21.1	15.0	18.6	23.3	8.5E-06	6.0E-02	0.0	9.7E-01	-2.8	5.9E-03
ENSG00000259139	DUX4L12**,***	0.6	1.2	0.7	22.7	1.2E-05	6.0E-02	-1.6	1.1E-01	2.7	8.0E-03
ENSG00000108691	CCL2ǂ	34.9	25.7	17.4	21.0	2.8E-05	6.5E-02	-4.1	3.8E-05	-5.4	7.2E-08
ENSG00000118503	TNFAIP3ǂ	76.1	56.4	41.3	20.9	3.0E-05	6.5E-02	-3.6	3.1E-04	-5.2	2.4E-07
ENSG00000162772	ATF3ǂ	113.3	79.6	54.3	21.6	2.0E-05	6.5E-02	-4.1	4.5E-05	-5.2	2.0E-07
ENSG00000169926	KLF13	1410.0	1957.8	1319.0	21.2	2.5E-05	6.5E-02	-1.3	2.0E-01	1.8	6.5E-02
ENSG00000179943	FIZ1	112.1	171.3	106.7	21.3	2.4E-05	6.5E-02	-2.2	2.8E-02	1.1	2.9E-01
ENSG00000204388	HSPA1Bǂ	5785.3	4663.8	2897.9	21.0	2.7E-05	6.5E-02	-4.4	1.2E-05	-4.8	2.0E-06
ENSG00000204389	HSPA1Aǂ	5660.5	4510.7	2819.8	20.7	3.2E-05	6.5E-02	-4.4	1.4E-05	-4.8	1.8E-06
ENSG00000206530	WDR52	333.8	307.9	345.8	20.6	3.3E-05	6.5E-02	1.8	7.4E-02	-1.0	3.3E-01
ENSG00000213058	RP4-765C7.2***	5.2	3.8	4.1	21.0	2.7E-05	6.5E-02	-0.7	4.6E-01	-3.3	1.1E-03
ENSG00000235950	DUX4L15**,***	1.2	2.5	1.3	20.6	3.4E-05	6.5E-02	-1.8	7.6E-02	2.6	9.3E-03
ENSG00000237919	RP11-181B18.1***	56.1	41.5	50.6	21.1	2.7E-05	6.5E-02	0.8	4.1E-01	-2.0	4.1E-02
ENSG00000259034	DUX4L3**	1.1	2.2	1.1	20.5	3.6E-05	6.5E-02	-2.0	4.3E-02	2.3	2.2E-02
ENSG00000259128	DUX4L2**	1.4	2.8	1.4	20.5	3.6E-05	6.5E-02	-2.2	2.8E-02	2.0	4.2E-02
ENSG00000006327	TNFRSF12Aǂ	10.8	8.7	6.5	19.9	4.8E-05	6.9E-02	-2.9	3.3E-03	-4.4	8.8E-06
ENSG00000125538	IL1Bǂ	37.5	25.6	20.1	19.9	4.8E-05	6.9E-02	-3.0	2.7E-03	-4.9	1.0E-06
ENSG00000170345	FOSǂ	805.4	562.7	384.1	19.9	4.8E-05	6.9E-02	-4.0	5.4E-05	-4.9	8.2E-07
ENSG00000198805	PNPǂ	78.7	57.9	46.3	20.0	4.4E-05	6.9E-02	-3.1	1.8E-03	-4.7	2.6E-06
ENSG00000259063	DUX4L5**	1.1	2.1	1.1	20.1	4.3E-05	6.9E-02	-2.0	4.5E-02	2.3	2.3E-02
ENSG00000124762	CDKN1A*,ǂ	134.3	112.1	77.3	19.5	5.8E-05	7.8E-02	-3.6	3.7E-04	-4.4	8.8E-06
ENSG00000146063	TRIM41	372.8	592.9	363.9	19.5	5.9E-05	7.8E-02	-1.6	1.1E-01	1.6	1.2E-01
ENSG00000236197	AC002429.5***,ǂ	51.1	37.9	45.8	19.3	6.5E-05	8.3E-02	0.3	8.0E-01	-2.4	1.7E-02
ENSG00000258834	DUX4L4*,***	1.1	2.2	1.2	19.2	6.8E-05	8.5E-02	-1.8	6.8E-02	2.4	1.7E-02
ENSG00000132002	DNAJB1ǂ	1361.4	1061.8	722.1	18.8	8.4E-05	9.2E-02	-3.9	8.7E-05	-4.8	2.0E-06
ENSG00000152580	IGSF10ǂ	66.6	52.6	63.3	18.8	8.4E-05	9.2E-02	1.5	1.3E-01	-1.2	2.3E-01
ENSG00000183018	SPNS2	348.0	536.9	366.0	18.9	7.8E-05	9.2E-02	-0.1	9.2E-01	2.9	3.8E-03
ENSG00000241449	RP11-545G3.1	22.7	17.0	20.9	18.8	8.3E-05	9.2E-02	0.6	5.4E-01	-2.1	3.9E-02

In the LRT, these genes survived corrected p<0.1 when applying the default DESeq2 Cooks Cutoff for outlier replacement (with trimmed mean), however they did not survive P<0.1 when a more conservative Cooks cutoff (1.25) was applied.

denotes GrCh37 Ensembl IDs which are deprecated (i.e. have been replaced by one or more new identifiers) in the current version (GrCh38).

denotes known unprocessed and/or antisense pseudogenes

denotes genes with adjusted p-values >0.1 FDR corrected when covarying for RIN scores in a reduced sample (see methods).

Figure 1

Top four differentially expressed genes among in MDD with and without suicide (Table 2, p<0.1 FDR corrected). Top left panel: Humanin-like 8 (MTRNR2L8) is higher in depression (MDD+MDD-S vs. CON, log2FC = 0.67, adjusted p=8.2E-5) and higher in suicide (MDD-S vs. CON+MDD, log2FC = 0.72, adjusted p=4.1E-9). Top right panel: Serpin peptidase inhibitor, clade H (heat shock protein 47), member 1 (SERPINH1) is lower in depression (log2FC=-0.72, adjusted p=2.5E-5) and lower in suicide (log2FC=-0.40). Lower left panel: Interleukin 8 (IL8) is lower in depression (log2FC=-0.53, adjusted p=0.0006). Lower right panel: Chemokine (C-C motif) ligand 4 (CCL4) is lower in depression (log2FC=-0.59, adjusted p=0.0003).

In the main LRT analysis, differential expression for one gene (ENSG00000088881, Gene name: early B-cell factor 4, Gene Symbol: EBF4) that did not survive p<0.1 adjusted became significant at p<0.1 adjusted in the analysis that covaried for RIN score (shaded row in Supplementary Table 3). This effect was driven by greater EBF4 expression in the MDD relative to both CON (p=1.2E-4 uncorrected) and MDD-S (p=1.27E-5) groups (data not shown). The top ten differentially expressed genes for the model that included RIN score, excluding deprecated GrCh37 Ensembl Ids and known unprocessed transcripts or pseudogenes, are listed in Supplementary Table 4.

The rightmost columns of Table 2 list the post hoc Wald statistics and uncorrected p-values for the depression and suicide contrasts (see methods). Of the differentially expressed genes listed in Table 2, all showed a depression-related effect (MTRNR2L8 was higher while the majority was lower in depression), while 7 of these genes also showed suicide effects which survived FDR < 0.1 (MTRNR2L8 was higher while the majority was lower in suicide). None of these genes appeared to show an effect specific to suicide (Table 2, second column from the right). The top gene (MTRNR2L8) showed higher expression in MDD (vs. CON, log2FC=0.46, p=0.02 uncorrected, data not shown). However, MTRNR2L8 showed highest mean (and variance) for the MDD-S group relative to both MDD (MDD-S vs. MDD, log2FC=0.87, p=8.2E-7 uncorrected, Supplementary Table 2) and CON groups (MDD-S vs. CON, log2FC=1.05, p=1.8E-10 uncorrected, Supplementary Table 2), suggesting that MDD-S (vs. MDD alone) represents a more severe and heterogeneous phenotype with respect to expression of this gene.

For miRNA, no differences survived the applied threshold for significance after correction for multiple comparisons (Kruskal-Wallis ANOVA with 3 groups: MDD-S, MDD and CON, p < 0.1 FDR corrected). Based on a recent review which identified specific miRNA species as playing a role in depression and/or suicide (38), we also restricted miRNA expression analyses to mir-185 and mir-491–3p using an uncorrected threshold, but did not observe group differences (p>0.2, data not shown).

Gene set analyses

Gene ontology (GO) functional enrichment analyses identified significant differences in 49 functional pathways at p<0.1 corrected (Table 3). The most significant GO gene sets include “chemokine receptor binding” (Figure 2), “cellular response to lipopolysaccharide” (Figure 3), and “positive regulation of angiogenesis, which were all lower in both depression and suicide comparisons groups, with suggestive evidence for lowest expression in the MDD-S group (Supplementary Table 2, see Supplementary Discussion). GO groups remained consistent when covarying for RIN scores (Supplementary Table 5).

Table 3

Significant gene ontology terms from ErmineJ using P-values for overall effect of condition (likelihood ratio test for ~ age + sex + condition vs. ~ age + sex, multifunction and FDR corrected p<0.1). ID=GO term numerical ID; Size=number of genes in GO term. Score=ErmineJ derived score for that GO term, Pval=uncorrected pvalue, Corr Pval=FDR corrected pvalue; MFCorrPvalue= multifunctionality FDR corrected pvalue. Correlation Pval=pvalue for the correlation amongst all genes in the GO term across all samples. For a detailed explanation of the rightmost five columns see (31,32).

Name	ID	Size	Score	Pval	Corr Pval	MFCorr Pvalue	Correlation Pval
chemokine receptor binding	GO:0042379	33	1.36	1.00E-12	3.32E-09	1.11E-09	1.00E+00
cellular response to lipopolysaccharide	GO:0071222	86	0.95	9.25E-11	2.04E-07	2.54E-06	1.00E+00
positive regulation of angiogenesis	GO:0045766	76	0.96	1.27E-10	2.10E-07	1.32E-03	1.00E+00
cellular response to biotic stimulus	GO:0071216	100	0.91	6.99E-10	9.26E-07	1.92E-05	1.00E+00
negative regulation of multi-organism process	GO:0043901	97	0.91	8.98E-10	9.92E-07	3.10E-03	1.00E+00
positive regulation of vasculature development	GO:1904018	84	0.92	1.42E-09	1.34E-06	1.17E-02	1.00E+00
cellular response to molecule of bacterial origin	GO:0071219	91	0.92	2.84E-09	2.36E-06	2.11E-05	1.00E+00
response to interleukin-1	GO:0070555	57	0.99	3.38E-09	2.49E-06	4.49E-04	1.00E+00
cellular response to interleukin-1	GO:0071347	37	1.07	2.19E-08	1.45E-05	1.70E-05	1.00E+00
positive regulation of chemotaxis	GO:0050921	78	0.90	3.86E-08	2.33E-05	2.93E-08	1.00E+00
serine-type endopeptidase inhibitor activity	GO:0004867	65	0.94	5.80E-08	3.21E-05	6.63E-09	1.00E+00
collagen metabolic process	GO:0032963	65	0.93	1.29E-07	6.57E-05	1.33E-09	1.00E+00
multicellular organismal macromolecule metabolic process	GO:0044259	70	0.90	1.67E-07	7.92E-05	8.29E-10	1.00E+00
regulation of cytokine biosynthetic process	GO:0042035	69	0.89	1.93E-07	8.53E-05	2.35E-02	1.00E+00
multicellular organismal metabolic process	GO:0044236	72	0.89	2.54E-07	1.05E-04	9.47E-10	1.00E+00
positive regulation of cytokine biosynthetic process	GO:0042108	43	0.99	5.31E-07	2.07E-04	2.16E-02	1.00E+00
positive regulation of leukocyte chemotaxis	GO:0002690	49	0.95	8.42E-07	3.10E-04	5.87E-08	1.00E+00
positive regulation of leukocyte migration	GO:0002687	66	0.90	9.81E-07	3.42E-04	6.05E-07	1.00E+00
negative regulation of extrinsic apoptotic signaling pathway	GO:2001237	69	0.87	1.21E-06	4.02E-04	1.08E-04	1.00E+00
regulation of tumor necrosis factor production	GO:0032680	62	0.89	1.64E-06	5.18E-04	7.55E-03	1.00E+00
regulation of leukocyte migration	GO:0002685	89	0.84	2.59E-06	7.16E-04	5.87E-05	1.00E+00
regulation of extrinsic apoptotic signaling pathway in absence of ligand	GO:2001239	46	0.93	2.89E-06	7.67E-04	2.32E-02	5.00E-04
regulation of tumor necrosis factor superfamily cytokine production	GO:1903555	64	0.87	4.38E-06	1.12E-03	1.38E-02	1.00E+00
positive regulation of behavior	GO:0048520	89	0.83	4.73E-06	1.16E-03	1.48E-05	1.00E+00
granulocyte chemotaxis	GO:0071621	34	1.00	5.09E-06	1.20E-03	1.57E-06	1.00E+00
regulation of leukocyte chemotaxis	GO:0002688	59	0.88	5.86E-06	1.34E-03	1.45E-06	1.00E+00
regulation of interleukin-2 production	GO:0032663	32	0.99	7.88E-06	1.69E-03	1.13E-02	1.00E+00
collagen binding	GO:0005518	49	0.91	9.03E-06	1.82E-03	3.32E-09	1.00E+00
response to temperature stimulus	GO:0009266	95	0.80	1.12E-05	1.95E-03	9.78E-02	1.00E-12
granulocyte migration	GO:0097530	35	0.98	1.11E-05	1.98E-03	4.37E-06	1.00E+00
leukocyte chemotaxis	GO:0030595	69	0.83	1.22E-05	2.03E-03	5.96E-05	1.00E+00
cellular response to tumor necrosis factor	GO:0071356	70	0.83	1.26E-05	2.04E-03	1.16E-03	1.00E+00
myeloid leukocyte migration	GO:0097529	52	0.89	1.48E-05	2.28E-03	7.45E-05	1.00E+00
response to tumor necrosis factor	GO:0034612	84	0.80	2.95E-05	3.83E-03	1.36E-02	1.00E+00
collagen fibril organization	GO:0030199	35	0.95	3.78E-05	4.64E-03	1.66E-09	9.60E-04
collagen biosynthetic process	GO:0032964	5	3.63	6.50E-05	7.07E-03	2.31E-02	8.31E-02
activation of signaling protein activity involved in unfolded protein response	GO:0006987	55	0.83	1.82E-04	1.53E-02	9.47E-05	7.00E-05
positive regulation of nuclease activity	GO:0032075	57	0.81	1.81E-04	1.54E-02	6.04E-05	7.00E-05
endoplasmic reticulum unfolded protein response	GO:0030968	73	0.78	3.63E-04	2.67E-02	2.37E-04	5.00E-06
cellular response to unfolded protein	GO:0034620	73	0.78	3.63E-04	2.67E-02	2.37E-04	5.00E-06
ER-nucleus signaling pathway	GO:0006984	79	0.76	3.87E-04	2.82E-02	4.41E-04	5.00E-06
regulation of nuclease activity	GO:0032069	61	0.80	4.46E-04	3.18E-02	3.56E-04	2.00E-05
skeletal muscle cell differentiation	GO:0035914	41	0.84	5.08E-04	3.44E-02	5.99E-02	1.00E+00
CCR chemokine receptor binding	GO:0048020	10	1.76	6.75E-04	4.34E-02	9.59E-02	1.00E+00
endoplasmic reticulum- Golgi intermediate compartment	GO:0005793	60	0.78	9.38E-04	5.50E-02	2.21E-09	7.00E-05
cellular response to topologically incorrect protein	GO:0035967	78	0.74	1.45E-03	6.78E-02	7.30E-04	5.00E-06
unfolded protein binding	GO:0051082	75	0.73	2.85E-03	9.66E-02	7.37E-10	1.00E+00
blood microparticle	GO:0072562	79	0.70	9.29E-03	1.93E-01	3.95E-02	1.00E+00

Figure 2

Details for GO gene set “chemokine receptor activity”

Each row corresponds to a gene that is included in the same set. “Element” column lists ensemble ID, “score” refers to the uncorrected p-value for the likelihood ratio test (LRT) when comparing the full model with condition (three levels: CON, MDD-S and MDD) vs. the reduced model (~ age + sex + condition vs. ~ age + sex, see methods), “QQ Score” graphs the scores: blue lines represent the observed scores, light grey line shows the expected distribution based on chance, “multifunctionality” indicates the multifunctionality of the gene. The number in parentheses is the number of annotations (e.g., GO terms) the gene has, which is roughly proportional to the multifunctionality, but the exact multifunctionality score takes into account the size of the groups to which the gene belongs, “QQ Multifunct” is similar to the QQ Score column but for the multifunctionality. If the gene set has “typical” multifunctionality, the pink line will tend to be near the grey line.

Figure 3

Details for GO gene set “cellular response to lipopolysaccharides”

See Figure 2 legend for description of columns.

In both suicide and depression contrasts (see methods), there was lower expression of transcripts involved in immune-related and microglial cellular functions (i.e. cellular response to type 1 interferon, regulation of leukocyte chemotaxis, regulation of cytokine biosynthetic process, etc. Table 3, Supplementary Table 6 and 7, bottom halves). None of the main results from Table 3 (LRT testing for overall effect of modeling group status) included GO sets which were higher in suicide or depression. Several GO terms appear to be more expressed higher in suicide and/or depression in exploratory contrast analyses: however, many of these results are driven by relatively weak effects, and all but a few consist of highly inter-correlated genes (Supplementary Tables 6 and 7, top halves).

Hypothesis-driven GO pathway analyses in suicide and depression

In addition to exploratory, whole-exome GO analyses above, we examined specific gene sets and pathways which were identified by published studies that used oligonucleotide microarrays to examine brain expression differences in MDD-S, CON and non-depressed suicides (8). In contrast to exploratory analyses, which were stringently corrected for multiple comparisons, we report uncorrected p-values for these hypothesis-based analyses.

Glial cell pathology

Previous studies implicate a central glial cell pathology and/or reduction in suicide and depression (8,39,40). Here, GO analyses were restricted to gene-sets which contained the keywords “oligodendrocyte” and “astrocyte”. Altered expression of genes with functions related to ‘oligodendrocyte differentiation’ (p=0.004, Supplementary Figure 1A) was observed, with lower expression in depression (p=0.016) but not suicide (p=0.17, data not shown). Altered expression of genes involved in ‘astrocyte cell migration’ (p=0.002, Supplementary Figure 1B) was observed, with lower expression in depression (p=1.5E-05) and in suicide (p=0.002, data not shown).

GABAergic and glutamatergic neurotransmission

Previous findings suggest altered expression of genes involved with GABAergic and glutamatergic neurotransmission in depression and suicide (8). The glutamatergic NMDA receptor antagonist ketamine appears to improve depressive symptoms much faster than traditional antidepressants, and greater expression of glutamatergic receptors in dlPFC has been observed in MDD-S (41). Here, we did not observe altered expression of transcripts in gene sets related to GABAergic neurotransmission (Supplementary Figure 1C), however we found less expression of transcripts involved in ‘regulation of synaptic transmission, glutamatergic’ (p=0.03, Supplementary Figure 1D), in both depression (p=0.05) and suicide (p=0.05, data not shown). The top gene in this GO set was the oxytocin receptor (OXTR) whose expression was lower in both depression (p=0.008) and suicide (p=0.001, data not shown).

ATP biosynthesis and ATPase coupled transmembrane activity

A previous study found altered expression of genes involved in adenosine 5′-triphosphate (ATP) production in the tricarboxylic acid (TCA) cycle and ATPase activity coupled to transport of ions across the cell membrane (8). Here we did not observe evidence for dysregulation in these pathways in suicide or depression (all p>0.7, Supplementary Figure 1E). However, we note that genes involved in ‘DNA-dependent ATPase activity’ were higher in the suicide comparison but not in the depression comparison (p<0.1 corrected, Supplementary Table 5). A post-hoc contrast in MDD-S vs. MDD groups only (no CON included) confirmed that this increase was specific to suicide (p=3.6E-06 uncorrected, Supplementary Figure 2).

Differential exon-usage in suicide and depression

We applied DEXSeq (19) to whole-exome RNA-seq data (see methods) in order to test for differential exon usage (DEU, i.e. differential splicing) in MDD with or without suicide. An exon in the gene for ATPase, class II, type 9B (ATP9B) was significant at p<0.1 after correction for multiple hypothesis testing (FDR adjusted p=0.087, Table 4, first row, Figure 5). Transcripts (splice variants) for the ATP9B gene in reference to exon bin 21 are shown in Supplementary Figure 3. As an exploratory analysis, we list the top 10 genes showing evidence for DEU in Table 4. These genes include microtubule associated serine/threonine kinase 3 (MAST3, adjusted p=0.14) and eukaryotic translation initiation factor 3, subunit K (EIF3K, adjusted p=0.36).

Table 4

Top 10 genes with evidence for differential exon usage in MDD and suicide

Mean normalized reads counts mapping to a particular exon bin (exon ID) for each group are listed in columns CON, MDD and MDD-S. LRT=Likelihood Ratio Test for comparing ~ sample + exon + exon:condition + age:exon + sex:exon as the full model vs. ~ sample + exon + age:exon + sex:exon as the reduced model (see DEXSEQ methods). Pvalue=uncorrected pvalue, padj=FDR corrected pvalue. Start/End=genomic locus start/end points in base pair units for a particular exon bin (exon ID). Width=exon length in base pairs.

Ensembl ID	GeneSymbol	Exon ID	BaseMean	dispersion	LRT(χ2)	pvalue	padj	CON	MDD	MDD-S	Start	End	Width
ENSG00000166377	ATP9B	E021	58.500	0.015	26.973	1.39E-06	0.083	12.531	14	13.232	76886267	76886375	109
ENSG00000099308	MAST3	E013	143.839	0.005	24.474	4.85E-06	0.144	18.350	17.030	18.229	18239210	18239268	59
ENSG00000178982	EIF3K	E024	164.541	0.001	21.859	1.79E-05	0.355	18.981	18.242	19.007	39125633	39125742	110
ENSG00000065609	SNAP91	E024	72.070	0.015	21.199	2.49E-05	0.371	NA	NA	NA	84302901	84302978	78
ENSG00000220997+ ENSG00000241621+ ENSG00000107951	GOLGA2P6, MTPAP	E009	59.977	0.002	20.615	3.34E-05	0.397	13.722	13.096	12.846	30625732	30625956	225
ENSG00000141504	SAT2	E020	59.199	0.001	19.664	5.37E-05	0.533	13.071	13.934	12.841	7530461	7530544	84
ENSG00000229896+ ENSG00000095951	HIVEP1	E026	105.764	0.004	18.899	7.87E-05	0.549	16.231	16.875	15.698	12163516	12164808	1293
ENSG00000116857	TMEM9	E005	95.882	0.001	18.822	8.18E-05	0.549	15.856	15.233	15.966	201104752	201104910	159
ENSG00000128335	APOL2	E007	55.346	0.003	18.792	8.31E-05	0.549	13.028	12.095	12.939	36629448	36629501	54
ENSG00000118473	SGIP1	E033	77.427	0.023	18.524	9.50E-05	0.565	NA	NA	NA	67147552	67147976	425

Figure 5

Exon-level expression in ATP9B after adjusting for group differences in gene-level expression

Top plot shows mean normalized counts for each of 50 splicing events (i.e. bins, or “exons”, see methods) in the ATP9B gene across for the three groups CON, MDD and MDD-S. There was a significant group effect in mean normalized counts in exon 21 (circled, adjusted p=0.087, see Table 1) after removing overall effect of gene expression and correcting for multiple comparisons across the whole exome. Blue lines indicate differentially expression exons surviving a more lenient threshold (adjusted p<0.1 correcting for only the 50 tests within this gene). The bottom row depicts 50 of 93 total unique splicing events (here referred to as “exon”) along the ATP9B gene. Exon bin 21 in reference to the 26 known ATP9B transcripts (isoforms) according to Human genome assembly GRCh37 is shown in Supplementary Figure 2.

Structural gene variant association analysis

An exploratory gene association analysis was conducted in search of novel or rare (primarily exonic) SNPs or insertions/deletions (indels) that may be of interest in terms of altered expression or damaging non-synonymous mutations (see RNA-seq Variant Discovery in Methods). We assessed case-control differences in allele frequencies in SUI vs NonSUI and MDD vs NonMDD groups. No gene variants survived correction for ultiple comparisons (FDR <0.1). In an exploratory analysis, variants in the top 34 genes listed in Table 2 reached threshold at p<0.05 uncorrected, while variants from whole-transcriptome reached threshold at p<0.0001 uncorrected. Supplementary Table 8 lists the putative effects and relative impact (i.e. non-synonymous mutation, exon-ablation etc. would have HIGH impact) of these variants according to snpEff v4.1 (36). None of these variants were predicted to have MODERATE or HIGH impact. Due to the limited number of samples with called variants (many samples were missing calls likely due in part to differences in sequencing depth across samples, see Supplementary Table 8) we did not attempt eQTL analyses.

Confirmation of differential expression results using available microarray data and qPCR

We sought further confirmation of RNAseq differential expression findings using a subset of subjects (2 CON and 3 MDD-S) for whom microarray data were also available from a previously published study from our group (26). This previous microarray study compared MDD-S vs. CON groups, thus replication analyses were restricted to those genes reported for MDD-S vs. CON differences listed in Supplementary Table 2. Briefly, probeset signal intensities were extracted with the robust multi-array average (RMA) algorithm (42) and averaged for each gene. There were 23 total genes identified by RNAseq (Supplementary Table 2), all having greater expression in CON vs. MDD-S, that were also included on the Affymetrix U133A platform. Due to limited sample size we did not attempt to conduct statistical inference within each gene. Instead, we tested for overall agreement in the directionality of findings across the 23 genes using a sign-test over their effect sizes for CON vs. MDD-S differences. Consistent with our RNA-seq findings, the majority of genes (17/23) had greater mean expression in the CON vs. MDD-S group (sign-test p=0.03, see Supplementary Table 9).

In addition, we sought further confirmation of results in the top 5 genes using qPCR over the full sample dataset. In general, qPCR results were consistent with RNA-seq results with the exception of MTRNR2L8 (see Supplementary Table 10).

Discussion

Exploratory whole-exome analysis revealed higher expression of humanin-like 8 (MTRNR2L8, or HN8) in both suicide and depression (Figure 1, Supplementary Table 2). MTRNR2L8 is an isoform of the humanin gene (HN) (43), a newly discovered mitochondrial-derived 24 amino acid peptide (44), which is coded by nuclear DNA. According to Uniprot, MTRNR2L8 plays a role as a neuroprotective and antiapoptotic factor. Accumulating evidence indicates HN plays a strong role in stress resistance through its neuroprotective, anti-apoptic and anti-inflammatory properties (45). Here, increased expression of MRTNR2L8 suggests a compensatory and constitutive up-regulation in response to the chronic stress in MDD.

With the exception of MRTNR2L8, differentially expressed genes (p<0.1 FDR corrected, Tables 3 and 4) had lower expression in depression and/or suicide (Supplementary Tables 5 and 6). This observation is consistent with reports of increased DNA methylation observed in depressed suicides (46). Studies that combine methylation and gene expression data will help further clarify the relationship between DNA methylation, gene expression and psychopathology (46).

Spermidine/spermine N-1 acetyltransferase (SAT1, or SSAT) is associated with suicide based on microarrays and has lower brain expression postmortem in depressed suicides (3,24,47), a finding which has been replicated in this dataset at more lenient threshold (p<0.005 uncorrected). See (48) for isoform-level profiling of the SAT1 gene and putative miRNA regulators (49) in the current sample. See Supplementary Material for further discussion of hypothesis-driven gene ontology and miRNA expression results.

Altered immune-related gene expression in depression and suicide

The most pronounced brain expression differences in depressed subjects with and without suicide were found in molecular pathways involved in microglial and immune system functions (in particular genes for chemokine receptors and ligands such as chemokine ligand 2 and 4 (CCL2 and CCL4) (50)) and expression was lower in MDD and MDD-S groups. Transcripts for Interleukin-8 (IL8, or CXCL8), a leukocyte chemotactic activating, pro-inflammatory cytokine (chemokine) (51), showed the greatest deficit in MDD relative to CON. No separate suicide effect was detected. Our findings of altered immune related gene expression are consistent with a closely related microarray gene expression study which found pathways involved in “immune cell activation” in BA46 (a region in frontal cortex adjacent to BA9) were altered in MDD and suicide (Table 1 in (8)).

Although cytokines do not appear to be essential mediators of depressive illness (52), they are thought to play a role in immune-mediated depressive disorders and may contribute to depressive symptoms in some patients (53). While many studies have found higher expression of pro-inflammatory cytokines such as IL6 in MDD (54), others have reported lower serum levels of IL6 (55). Other studies have found that plasma and CSF levels of IL-8 and chemokines were lower in suicide attempters compared to non-attempters (56), while one study found that proinflammatory cytokines such as IL6 and TNF-a were higher in the brain (BA10) of teenage suicides(57)

Others have found less expression of immune-related genes such as CD44, MARCH1 and CD300LB in frontal cortex of depressed suicides (58). A recent study used RNA-seq to examine whole-exome expression differences in blood from a large live sample and found higher expression of type I interferon signaling in MDD compared with healthy volunteers (59). Our finding of lower chemokine and immune related gene expression in MDD-S appear to contradict recent findings of increased microglia priming and chemokine gene expression in white matter (dorsal anterior cingulate) of MDD-S postmortem (60) although are consistent with recent meta-analyses and studies which found lower levels of pro-inflammatory cytokines such as IL8 (61–63) as well as chemokines such as CCL2 and CCL4 (64) in CSF of suicidal patients vs. non-suicidal patients and controls. Despite robust experimental evidence indicating that elevated chemokine expression and subsequent recruitment of macrophages into the brain is detrimental, increasing evidence suggests that chemokines also possess pleiotropic and beneficial properties beyond chemotaxis (50). CCL2 may be important in neurogenesis as recombinant CCL2 was first shown to promote glial cell proliferation and growth in vitro (65), while recent evidence suggests that CCL2 may direct differentiation of precursor cells into neurons, astrocytes and oligodendrocytes (66). Additional neuroprotective and neurotrophic functions of the ligands CXCL8, CXCL1, and CXCL2 have been shown, which include mediating self-defense mechanisms against Fas-initiated apoptotic cell death, and modulating synaptic transmission through altering calcium channel excitability and neurotransmitter release (50). Other lines of research show that microglia serve important functions in the uninjured brain through pruning of unwanted synapses, (67) and it has been shown that primary deficits in microglia lead to deficient synaptic pruning, decreased functional brain connectivity, deficits in social interaction and increased repetitive behavioral phenotypes in mice (68).

Our results show marked lower expression of genes associated with microglia cell functions (chemokine activity, regulation of cytokine biosynthetic process) in depression, suggesting that altered microglia-mediated synaptic pruning may be present in MDD (68). In addition, mounting evidence suggest reductions in astrocytes and aberrant astrocyte function in MDD and MDD-S (69–71). Since chemokine and immune related genes are expressed in microglia and astrocytes during normal physiological conditions and play a role in neuron-glia communication (72,73) and neuroprotection (74), our findings of reduced chemokine and immune-related gene expression may reflect cellular abnormalities (i.e. reduced astrocytes) present in the brain of MDD-S. Indeed, three of the top 4 genes listed in Table 1 appear to be expressed specifically in microglia and/or astrocytes according to gene expression data from mouse cortex (http://oganm.github.io/celltypes, Ogan Mancarci, Lilah Toker and Paul Pavlidis, personal communication, see Supplementary Figure 4). Moreover, our findings are consistent with a recent PET study examining neuroinflammation marker translocator protein (TPSO) and which failed to find evidence for increased neuroinflammation in mild to moderate depressed individuals (in fact 70% of the depression sample showed lower [11C]PBR28 binding relative to genotype matched controls (75). PET brain imaging of microglial activation after lipopolysaccharide administration (76) in larger samples of depressed and/or suicidal individuals may help clarify whether increased neuroinflammatory activity is associated with these phenotypes. Lower immune-related gene expression in MDD may also depend on tissue-type (brain vs. blood) as well as location in the brain (77). Future studies are required to further understand the function of pro-inflammatory cytokines and immune response pathways in cerebral cortex and their role in depression and suicide.

Reduced angiogenesis gene expression in depression and suicide

Chemokines are also known to play an important role in the regulation of angiogenesis (78). Gene ontology analysis revealed lower expression of other genes involved in regulation of angiogenesis (Supplementary Table 5, Figure 4). Although little is known about the role of angiogenesis in depression, hippocampal angiogenesis has been shown to be coupled to neurogenesis, and both are greater in antidepressant-treated major depression (79). Our results are consistent with some previous studies which showed high serum levels of endostatin, an angiogenesis inhibitor, were associated with late-life depression (80). These findings warrant future studies to further understand the role of cortical angiogenesis in the dorsal prefrontal cortex in major depression.

Figure 4

Details for GO gene set “positive regulation of angiogenesis”

See Figure 2 legend for description of columns.

Here we observed lower expression as well as evidence for alternative splicing of SERPINH1 in MDD and MDD-S (Figure 1, Table 3). SERPINH1, a member of the serine proteinase inhibitor gene family, encodes HSP47, which was initially identified as a gelatin binding protein named colligin (81). Functional polymorphisms and/or missense mutations in the SERPINH1 gene have been linked to preterm premature ruptures of membranes in African Americans (82) and osteogenic disorders (83). Although a putative role for this gene in depression and/or suicide remains unclear, GO analysis revealed lower expression of other genes involved in ‘serine-type endopeptidase inhibitor activity’, in both MDD and MDD-S groups (Table 4, Supplementary Tables 5 and 6). Serine endopeptidases are ubiquitous enzymes that cleave peptide bonds in serine proteases, in which serine serves as the nucleophilic amino acid at the enzyme's active site (84). Serine proteases are known to play a regulatory role in angiogenesis and vasculogenesis by activation, liberation and modification of angiogenic growth factors, degradation of the endothelial and interstitial matrix, and modification of the properties of angiogenic growth factors and cytokines (85). Excessive proteolysis causes loss of endothelial cell-matrix interaction and impairs angiogenesis. Thus lower expression and/or alternative splicing of serine peptidase inhibitor activity could potentially account for reduced expression of angiogenic growth factors observed here in the MDD and MDD-S groups.

Differential exon usage analyses

Interestingly, the gene with the strongest evidence for differential splicing (alternate exon-usage) in MDD and suicide was ATPase, class II type 9B (p<0.1 corrected, Table 4, Figure 5). In particular, exon bin 21 (corresponding to exon 5 of the primary protein coding transcript ATP9B-001, see Supplementary Figure 2) was higher in both MDD and MDD-S groups, suggesting a relationship to MDD, leaving unclear whether there is any specific relationship to suicide. This gene belongs to several GO categories previously implicated in depressed suicides, which include “ATPase activity, coupled to transmembrane movement of substances” and “hydrolase activity, acting on acid anhydrides, catalyzing transmembrane movement of substances” (8). ATP9B is a Type IV P-type ATPases (P4-ATPase), a class of putative phospholipid flippases that translocate phospholipids from the exoplasmic (lumenal) to the cytoplasmic leaflet of lipid bilayers (86). P4-ATPases play diverse and essential roles in membrane and transport vesicle biogenesis, and deficiency in P4-ATPases has been linked to liver disease in humans (87). In addition to their putative role in establishing and maintaining plasma membrane asymmetry, P4-ATPases are thought to be involved in vesicle-mediated protein transport in the exocytic and endocytic pathways as well as non-vesicular intracellular trafficking of sterols (87). Relatively little is known about the functions of ATP9B, other P4-ATPases and their roles in disease (88), and more research is necessary to understand the pathophysiological consequences of alterations in and alternative splicing of the ATP9B gene and its putative role in MDD and suicide.

Limitations

The main limitation of this study is modest sample sizes, but it is difficult to obtain larger sample sizes of cases and controls assessed in such detail clinically, neuropathologically and known to be medication and drug-free. Additionally, we examined a group of MDD subjects who did not die by suicide and were not treated with antidepressants, a group exceedingly difficult to obtain that has not been previously examined in the literature. In this study we did not observed significant effects specific to suicide and not depression (i.e. MDD-S vs. MDD). This is most likely due in part to the small sample size of the MDD non-suicide group (N=9). However, for several of the top genes and pathways it appeared that expression mean difference (and variance) were higher in the MDD-S relative to the MDD group (i.e. MTRNR2L8 was higher and angiogenesis and chemokine receptor binding were lower in MDD-S vs. MDD, see Supplemental Discussion). We anticipate that larger sample sizes of well-matched psychiatric control samples (i.e. medication-free MDD without suicide, a hard to obtain group) in future studies should help parse out effects of suicide vs. depression.

A second limitation is that we only examined one brain region and did not distinguish between neurons and glia. It is difficult to obtain sufficient transcript for quantification from each cell type. qPCR validations of the top 5 genes were consistent with RNA-seq results with the exception of one gene (MTRNR2L8). Because qPCR uses an exponential amplification system that introduces different types of biases and variability that are not observed with direct RNA-seq, it is arguable that RNA-seq should be the gold standard for differential expression microarray and qPCR analyses (89,90). Although RNA-seq has been shown to be highly accurate and has been well validated by others (17,18,88) new findings from this study (for MTRNRL28 in particular, and new findings from any other single study) should be regarded as provisional until replicated in independent samples.

In the current study, gene expression was quantified using most or all known RNA that is transcribed from a given gene, and that includes non-coding and alternate protein-coding isoforms, while gene ontology analyses ascribed functions to each gene based on only one or few protein-coding products each gene. However, it should be noted that RNA products from the same gene could have different, or even opposing functions (14). Our knowledge of non-coding (ncRNA) and isoform-level protein function is currently very limited, but as this knowledge grows so will the ability to conduct functional pathway analyses at the isoform level (14). Furthermore, alterations in gene expression in the cerebral cortex of MDD cases who died by suicide or other means does not necessarily confirm a causal relationship between altered expression and MDD or suicidal behavior.

Finally, lack of replicability in case-control studies is often attributed to the phenotypic heterogeneity of both MDD and suicide (i.e. two individuals can both have a diagnosis of MDD but share only a few symptoms) which has contributed to lack of robustly replicated genetic loci contributing to the MDD (89,90) and has helped prompt the Research Domain Criteria (RDoC) initiative to study dimensions of behavior personality that cut across binary disorder labels (91)).

Conclusion

In the current study, whole-exome sequencing (RNA-seq) of dorsolateral prefrontal cortex (BA9) was used to examine group differences in brain gene and miRNA expression between depressed suicides, depressed non-suicides, and sudden death non-psychiatric controls. Exploratory and hypothesis-driven GO functional pathway analyses revealed lower (higher) expression of genes involved in immune response, microglial and glial cell functions, angiogenesis, regulation of glutamatergic neurotransmission, (anti-apoptosis and neuroprotection) in suicides with or without MDD, while differential exon usage analysis suggests alternative splicing in ATP9B, a P4-ATPase transmembrane flippase involved in vesicular membrane biogenesis and protein transport. These effects appear to be driven primarily by depression, with the exception of genes involved in ‘DNA-dependent ATPase activity’, which appear to be greater in suicide but not depression. Results suggest cortical gene expression related to microglial, endothelial and glial cell functions may contribute toward or be altered in MDD and suicide, and identify the putative genes and pathways which may play a functional role in these disorders.