Advertisement

Replicable and Coupled Changes in Innate and Adaptive Immune Gene Expression in Two Case-Control Studies of Blood Microarrays in Major Depressive Disorder

Open AccessPublished:July 06, 2017DOI:https://doi.org/10.1016/j.biopsych.2017.01.021

      Abstract

      Background

      Peripheral inflammation is often associated with major depressive disorder (MDD), and immunological biomarkers of depression remain a focus of investigation.

      Methods

      We used microarray data on whole blood from two independent case-control studies of MDD: the GlaxoSmithKline–High-Throughput Disease-specific target Identification Program [GSK-HiTDiP] study (113 patients and 57 healthy control subjects) and the Janssen–Brain Resource Company study (94 patients and 100 control subjects). Genome-wide differential gene expression analysis (18,863 probes) resulted in a p value for each gene in each study. A Bayesian method identified the largest p-value threshold (q = .025) associated with twice the number of genes differentially expressed in both studies compared with the number of coincidental case-control differences expected by chance.

      Results

      A total of 165 genes were differentially expressed in both studies with concordant direction of fold change. The 90 genes overexpressed (or UP genes) in MDD were significantly enriched for immune response to infection, were concentrated in a module of the gene coexpression network associated with innate immunity, and included clusters of genes with correlated expression in monocytes, monocyte-derived dendritic cells, and neutrophils. In contrast, the 75 genes underexpressed (or DOWN genes) in MDD were associated with the adaptive immune response and included clusters of genes with correlated expression in T cells, natural killer cells, and erythroblasts. Consistently, the MDD patients with overexpression of UP genes also had underexpression of DOWN genes (correlation > .70 in both studies).

      Conclusions

      MDD was replicably associated with proinflammatory activation of the peripheral innate immune system, coupled with relative inactivation of the adaptive immune system, indicating the potential of transcriptional biomarkers for immunological stratification of patients with depression.

      Keywords

      Depression and inflammation are often associated with one another. Depressive symptoms in a large population sample were significantly related to blood concentrations of C-reactive protein (CRP; odds ratio ∼1.8 for depressive symptoms in people with CRP > 3 mg/L vs. CRP < 1 mg/L) (
      • Wium-Andersen M.K.
      • Ørsted D.D.
      • Nielsen S.F.
      • Nordestgaard B.G.
      Elevated C-reactive protein levels, psychological distress, and depression in 73,131 individuals.
      ). Multiple case-control studies of major depressive disorder (MDD) have reported increased peripheral blood concentrations of CRP (Cohen’s d ∼0.50) and proinflammatory cytokines such as interleukin 6 (d ∼0.50) and tumor necrosis factor (d ∼0.40) in MDD (
      • Valkanova V.
      • Ebmeier K.P.
      • Allan C.L.
      CRP, IL-6 and depression: A systematic review and meta-analysis of longitudinal studies.
      ,
      • Haapakoski R.
      • Mathieu J.
      • Ebmeier K.P.
      • Alenius H.
      • Kivimaki M.
      Cumulative meta-analysis of interleukins 6 and 1β, tumour necrosis factor α and C-reactive protein in patients with major depressive disorder.
      ). The prevalence of comorbid depression is increased in many nonpsychiatric inflammatory disorders (
      • Matcham F.
      • Rayner L.
      • Steer S.
      • Hotopf M.
      The prevalence of depression in rheumatoid arthritis: A systematic review and meta-analysis.
      ).
      There is growing evidence for a causal effect of inflammation on depression. Peripheral inflammation precedes the emergence of depressive symptoms in longitudinal epidemiological studies (
      • Khandaker G.M.
      • Pearson R.M.
      • Zammit S.
      • Lewis G.
      • Jones P.B.
      Association of serum interleukin 6 and C-reactive protein in childhood with depression and psychosis in young adult life: A population-based longitudinal study.
      ) and in about 30% of patients receiving proinflammatory interferon-α treatment for hepatitis C (
      • Bull S.J.
      • Huezo-Diaz P.
      • Binder E.B.
      • Cubells J.F.
      • Ranjith G.
      • Maddock C.
      • et al.
      Functional polymorphisms in the interleukin-6 and serotonin transporter genes, and depression and fatigue induced by interferon-α and ribavirin treatment.
      ,
      • Hepgul N.
      • Cattaneo A.
      • Agarwal K.
      • Baraldi S.
      • Borsini A.
      • Bufalino C.
      • et al.
      Transcriptomics in interferon-α-treated patients identifies inflammation-, neuroplasticity- and oxidative stress-related signatures as predictors and correlates of depression.
      ). Experimental challenge with peripheral proinflammatory stimuli in animals robustly induces a syndrome of illness behavior and anhedonia that approximates depressive symptoms (
      • Dantzer R.
      • O’Connor J.
      • Freund G.
      • Johnson R.
      • Kelley K.
      From inflammation to sickness and depression: When the immune system subjugates the brain.
      ). Peripheral immune cells and cytokines are known to mediate signals across the blood-brain barrier by several mechanisms (
      • Banks W.A.
      • Kastin A.J.
      • Gutierrez E.G.
      Penetration of interleukin-6 across the murine blood-brain barrier.
      ). Activation of microglia can locally amplify the effects of even a weak peripheral proinflammatory signal on neuronal function and behavior (
      • Perry V.H.
      • Cunningham C.
      • Holmes C.
      Systemic infections and inflammation affect chronic neurodegeneration.
      ).
      These observations suggest that pharmacological disruption of peripheral proinflammatory signals could be therapeutically effective, at least for a subgroup of patients with depression. However, it is most unlikely that any single anti-inflammatory drug will prove to be superior to all existing treatments for all patients (
      • Miller A.H.
      • Raison C.L.
      Are anti-inflammatory therapies viable treatments for psychiatric disorders? Where the rubber meets the road.
      ). Only about a third of patients with depression have biological evidence for peripheral inflammation [e.g., CRP > 3 mg/L (
      • Raison C.
      • Capuron L.
      • Miller A.
      Cytokines sing the blues: Inflammation and the pathogenesis of depression.
      )], and an anti-inflammatory drug seems likely a priori to be most effective for an inflamed subgroup of patients with depression. There are many markers of the peripheral immune system that can be conveniently measured in venous blood samples from patients with depression, including cytokines, CRP, and other proteins; cell counts from flow cytometry; and gene transcription (
      • Lin E.
      • Tsai S.J.
      Genome-wide microarray analysis of gene expression profiling in major depression and antidepressant therapy.
      ).
      Transcriptional (messenger RNA [mRNA]) biomarkers have the potential advantages of assay stability (compared with cytokines) and target specificity (compared with cell counts). However, previous case-control studies of peripheral blood gene expression and MDD have been inconsistent (
      • Jansen R.
      • Penninx B.W.J.H.
      • Madar V.
      • Xia K.
      • Milaneschi Y.
      • Hottenga J.J.
      • et al.
      Gene expression in major depressive disorder.
      ,
      • Garbett K.A.
      • Vereczkei A.
      • Kálmán S.
      • Brown J.A.
      • Taylor W.D.
      • Faludi G.
      • et al.
      Coordinated messenger RNA/microRNA changes in fibroblasts of patients with major depression.
      ,
      • Mostafavi S.
      • Battle A.
      • Zhu X.
      • Potash J.B.
      • Weissman M.M.
      • Shi J.
      • et al.
      Type I interferon signaling genes in recurrent major depression: Increased expression detected by whole-blood RNA sequencing.
      ,
      • Glahn D.C.
      • Curran J.E.
      • Winkler A.M.
      • Carless M.A.
      • Kent J.W.
      • Charlesworth J.C.
      • et al.
      High dimensional endophenotype ranking in the search for major depression risk genes.
      ,
      • Shelton R.C.
      • Claiborne J.
      • Sidoryk-Wegrzynowicz M.
      • Reddy R.
      • Aschner M.
      • Lewis D.A.
      • et al.
      Altered expression of genes involved in inflammation and apoptosis in frontal cortex in major depression.
      ) (Supplemental Table S1). To date, the Netherlands Study of Depression and Anxiety (NESDA) (
      • Jansen R.
      • Penninx B.W.J.H.
      • Madar V.
      • Xia K.
      • Milaneschi Y.
      • Hottenga J.J.
      • et al.
      Gene expression in major depressive disorder.
      ) is the largest single case-control study of the peripheral blood transcriptome in MDD (882 current patients with MDD, 635 remitted patients with MDD, and 331 control subjects). The NESDA reported statistically significant (false discovery rate [FDR] = 10%) differential expression of 129 genes enriched for interleukin 6 and natural killer (NK) cell signaling pathways (
      • Blume J.
      • Douglas S.D.
      • Evans D.L.
      Immune suppression and immune activation in depression.
      ).
      We were primarily motivated by the hypothesis that MDD is associated with peripheral blood transcriptional markers of innate immune system activation (
      • Dantzer R.
      • O’Connor J.
      • Freund G.
      • Johnson R.
      • Kelley K.
      From inflammation to sickness and depression: When the immune system subjugates the brain.
      ,
      • Miller A.
      • Maletic V.
      • Raison C.
      Inflammation and its discontents: The role of cytokines in the pathophysiology of major depression.
      ). We were also concerned to focus on results that were more likely to replicate across case-control studies of gene expression in MDD. We report Affymetrix microarray data on 18,863 probes from two independently designed and conducted case-control studies of MDD. We used a Bayesian method to identify genes that were differentially expressed in both studies. Focusing on a consensus set of 165 genes, we investigated the functional significance of the genes that were differentially (over- or under-) expressed in cases compared with controls. We also explored the secondary hypothesis that innate immune system activation is coupled to relative inactivation of the adaptive immune system in patients with MDD (
      • Grosse L.
      • Hoogenboezem T.
      • Ambrée O.
      • Bellingrath S.
      • Jörgens S.
      • de Wit H.J.
      • et al.
      Deficiencies of the T and natural killer cell system in major depressive disorder: T regulatory cell defects are associated with inflammatory monocyte activation.
      ,
      • Snijders G.
      • Schiweck C.
      • Mesman E.
      • Grosse L.
      • de Wit H.J.
      • Nolen W.A.
      • et al.
      A dynamic course of T cell defects in individuals at risk for mood disorders.
      ,
      • Grosse L.
      • Carvalho L.A.
      • Wijkhuijs A.M.
      • Bellingrath S.
      • Ruland T.
      • Ambree O.
      • et al.
      Clinical characteristics of inflammation-associated depression: Monocyte gene expression is age-related in major depressive disorder.
      ).

      Methods and Materials

      Samples

      We analyzed data from two case-control studies of depression: the GlaxoSmithKline–High-Throughput Disease-specific target Identification Program (GSK–HiTDiP) study and the Janssen–Brain Resource Company (Janssen–BRC) study. Other aspects of these studies have been previously reported (
      • Domenici E.
      • Willé D.R.
      • Tozzi F.
      • Prokopenko I.
      • Miller S.
      • McKeown A.
      • et al.
      Plasma protein biomarkers for depression and schizophrenia by multi analyte profiling of case-control collections.
      ,
      • Muglia P.
      • Tozzi F.
      • Galwey N.W.
      • Francks C.
      • Upmanyu R.
      • Kong X.Q.
      • et al.
      Genome-wide association study of recurrent major depressive disorder in two European case-control cohorts.
      ,
      • Tilahun A.
      • Lin D.
      • Shkedy Z.
      • Geys H.
      • Alonso A.
      • Peeters P.
      • et al.
      Genomic biomarkers for depression: Feature-specific and joint biomarkers.
      ,
      • Zoon H.
      • Veth C.
      • Arns M.
      • Drinkenburg W.
      • Talloen W.
      • Peeters P.
      • et al.
      EEG alpha power as an intermediate measure between brain-derived neurotrophic factor Val66Met and depression severity in patients with major depressive disorder.
      ,
      • Liu Y.
      • Yieh L.
      • Yang T.
      • Drinkenburg W.
      • Peeters P.
      • Steckler T.
      • et al.
      Metabolomic biosignature differentiates melancholic depressive patients from healthy controls.
      ); demographic and clinical details on the samples are provided in Supplemental Table S2.

      GSK–HiTDiP

      This study was designed primarily to test for an association between genetic (DNA) variation and diagnosis of depression. Minimal sociodemographic and clinical data were collected, but microarray data were available for analysis (after quality control) from whole-blood samples stored for less than 6 years on a sample comprising 113 patients with MDD prospectively balanced for comorbid anxiety disorder [57 with generalized anxiety disorder diagnosed by the Mini-International Neuropsychiatric Interview (
      • Sheehan D.V.
      • Lecrubier Y.
      • Sheehan K.H.
      • Amorim P.
      • Janavs J.
      • Weiller E.
      • et al.
      The Mini-International Neuropsychiatric Interview (MINI): The development and validation of a structured diagnostic psychiatric interview for DSM-IV and ICD-10.
      ) and 56 without anxiety disorder] and 57 healthy control subjects. All participants provided informed consent in writing. The study was approved by an independent ethics review board.

      Janssen–BRC

      This study was designed primarily for biomarker discovery. Microarray data were available for analysis (after quality control) from whole-blood samples stored for less than 1 year on a sample comprising 94 patients with MDD (40 with generalized anxiety disorder diagnosed post hoc by the Mini-International Neuropsychiatric Interview and 54 without anxiety disorder) and 100 healthy control subjects. Additional data on melancholic symptom severity, anxiety, substance use, and body mass index (BMI) were available for patients with MDD. All participants provided informed consent in writing. The study was approved by an independent ethics review board.
      Whole-blood samples from both studies were analyzed using the Affymetrix Human Genome U133 plus 2.0 array. We applied identical quality control, normalizing, and annotation algorithms to the two datasets, resulting in the estimation of mRNA expression at each of 18,863 unique probes for each participant.

      Differential Expression Analysis

      To determine differential gene expression between cases and controls, we adopted the same strategy for both studies; see Supplemental Figure S1 for a schematic overview of the data analysis strategy. For each gene i = 1, …, 18,863, we fitted a linear regression model that included group (coded Gr; two-level factor, case/control), batch (B; two-level factor, 1/2), gender (Ge; two-level factor, male/female), age (Ag; continuous), and presence/absence of anxiety (An; two-level factor, 0/1) as covariates. Denoting samples by j = 1, …, n and the design matrix by X, the model is
      yij=βi0+k=12βikGrxjkGr+l=12βilBxjlB+m=12βimGexjmGe+βiAgxjAg+r=12βirAnxjrAn+εi,
      (1)


      where εiN(0, σi2). For identifiability, we imposed the contrasts βi1Gr=βi1B=βi1Ge=βi1An=0 on model parameters. The model was fitted using the R package limma (
      • Ritchie M.E.
      • Phipson B.
      • Wu D.
      • Hu Y.
      • Law C.W.
      • Shi W.
      • et al.
      limma powers differential expression analyses for RNA-sequencing and microarray studies.
      ). Subsequently, we tested the null hypothesis Hi0: βi2Gr = 0; that is, there is no difference in expression of the ith gene between the two groups using the moderated t statistic (
      • Smyth G.K.
      Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.
      ). For each of the 18,863 probes in each study, the p value was generated using the asymptotic approximate distribution (and also nonparametrically by a permutation algorithm; see Supplemental Tables S3 and S4).

      Combining p Values for Differential Gene Expression From Two Studies

      To identify MDD-related genes that were replicated in both the GSK–HiTDiP and Janssen–BRC datasets, we set the p-value threshold for significance of differential expression of each gene in each study to optimize in some sense the number of genes that were differentially expressed in both studies compared with the number of coincidental differences expected by chance. To do this, we used the method of Blangiardo and Richardson (
      • Blangiardo M.
      • Richardson S.
      Statistical tools for synthesizing lists of differentially expressed features in related experiments.
      ), implemented in the R package sdef (
      • Blangiardo M.
      • Cassese A.
      • Richardson S.
      sdef: An R package to synthesize lists of significant features in related experiments.
      ), and specified the p-value threshold as q2, which represents the largest (most lenient) p value < .05 for which there are at least twice as many significant case-control differences in common between the two studies as expected by chance.

      Gene Coexpression Network Analysis

      We used weighted gene correlation network analysis (
      • Zhang B.
      • Horvath S.
      A general framework for weighted gene co-expression network analysis.
      ) to construct a normative coexpression network in which nodes represent genes and weighted edges represent correlations between the expression of pairs of genes in healthy control subjects (
      • Zhang B.
      • Horvath S.
      A general framework for weighted gene co-expression network analysis.
      ). To maximize the amount of data available for this estimation, and to ensure that the resulting network was representative of both studies, we included the healthy control data from both the GSK–HiTDiP and Janssen–BRC studies (N = 157 in total). We used the consensus weighted gene correlation network analysis method to construct a weighted undirected graph that could be decomposed into a set of modules of strongly coexpressed genes (
      • Zhang B.
      • Horvath S.
      A general framework for weighted gene co-expression network analysis.
      ). Eigengenes were used to test for case-control differences in the gene expression profile on average within each module of the normative weighted gene correlation network analysis transcriptome (Supplemental Table S11).

      Gene Ontology Enrichment Analysis

      We performed enrichment analysis for various gene lists such as the list of genes differentially expressed in both studies and the list of genes comprising each module of the coexpression network. We used the R package topGO (

      Alexa A, Rahnenfuhrer J (2010): topGO: Enrichment analysis for gene ontology. R package, version 2.0.

      ) to determine whether these lists were significantly enriched for specific Gene Ontology (GO) terms (
      • Ashburner M.
      • Ball C.A.
      • Blake J.A.
      • Botstein D.
      • Butler H.
      • Cherry J.M.
      • et al.
      Gene Ontology: Tool for the unification of biology.
      ) using a stringent Bonferroni correction (q < .05) for multiple comparisons across all 10,124 terms, resulting in a p-value threshold of 4.94 × 10−6.

      Protein–Protein Interaction Mapping

      We used the Search Tool for the Retrieval of Interacting Genes/Proteins [http://string-db.org (
      • Szklarczyk D.
      • Franceschini A.
      • Wyder S.
      • Forslund K.
      • Heller D.
      • Huerta-Cepas J.
      • et al.
      STRING v10: Protein–protein interaction networks, integrated over the tree of life.
      )] to determine the network of known protein–protein interactions among the genes differentially expressed in both studies (
      • Franceschini A.
      • Szklarczyk D.
      • Frankild S.
      • Kuhn M.
      • Simonovic M.
      • Roth A.
      • et al.
      STRING v9: 1. Protein–protein interaction networks, with increased coverage and integration.
      ).

      Cell Type–Specific Gene Expression Analysis

      We investigated the genes differentially expressed in both studies in relation to a third independent microarray (Affymetrix) dataset designed to examine cell type–specific expression patterns (
      • Mabbott N.A.
      • Baillie J.K.
      • Brown H.
      • Freeman T.C.
      • Hume D.A.
      An expression atlas of human primary cells: Inference of gene function from coexpression networks.
      ). In particular, we focused on the expression of MDD-related genes across 37 sorted cell samples of the following immune cell types: erythroblasts (8 samples), monocytes (6 samples), monocyte-derived dendritic cells (5 samples), neutrophils (3 samples), B cells (4 samples), CD4+ T cells (5 samples), CD8+ T cells (5 samples), and NK cells (6 samples). We used the BioLayout Express3D software (
      • Freeman T.C.
      • Goldovsky L.
      • Brosch M.
      • Van Dongen S.
      • Mazière P.
      • Grocock R.J.
      • et al.
      Construction, visualisation, and clustering of transcription networks from microarray expression data.
      ,
      • Theocharidis A.
      • van Dongen S.
      • Enright A.
      • Freeman T.
      Network visualization and analysis of gene expression data using BioLayout Express3D.
      ) to visualize clusters or subgroups of genes that shared similar patterns of expression across different cell types.

      Results

      Combining Differential Expression Across Two Case-Control Studies of MDD

      Separate differential expression analyses of the GSK–HiTDiP and Janssen–BRC datasets yielded two lists of p values for the same set of 18,863 gene expression probes (Supplemental Tables S3 and S4). The Bayesian method identified q2 = .0246 as the largest p value threshold associated with twice the number of genes differentially expressed in both studies compared with the number of coincidental case-control differences expected by chance. At this threshold, 173 genes were differentially expressed in both studies. We further refined this gene list by restricting attention to the 165 consensus genes (95%) that showed the same direction of overexpression (UP) or underexpression (DOWN) in the two case-control studies (Supplemental Table S5).
      These findings were corroborated by Fisher’s chi-square test for combining p values, which identified 393 genes that were differentially expressed in both studies with FDR = 10%, of which 146 were included in the list of replicably, concordantly, differentially expressed genes defined by Bayesian analysis (Supplemental Table S6).
      Furthermore, the set of 165 replicably concordant genes (henceforth the MDD-165 consensus set) was partially corroborated by the prior results of the NESDA case-control study (
      • Jansen R.
      • Penninx B.W.J.H.
      • Madar V.
      • Xia K.
      • Milaneschi Y.
      • Hottenga J.J.
      • et al.
      Gene expression in major depressive disorder.
      ). Considering the 15,830 genes that were measured by the microarrays used in all three studies (NESDA, Janssen–BRC, and GSK–HiTDiP), 150 of the MDD-165 consensus set were included, of which 7 genes (CD247, PRKCH, PGLYRP1, NFATC2, ST6GAL1, MAPK14, and MTSS1) were also differentially expressed, with the same sign of fold change, in the NESDA study at FDR = 5%. When the FDR threshold for the NESDA study was relaxed to 10% and 20%, 18 and 45 of the 150 genes differentially expressed in both the GSK–HiTDiP and Janssen–BRC studies, respectively, were also differentially expressed in the NESDA study (Supplemental Table S5). By comparison, we expect by chance to find 1.95, 8.30, and 30.40 genes in common with the NESDA study at FDR thresholds of 5%, 10%, and 20%, respectively. In contrast, there was less overlap between the MDD-165 consensus set and the top 29 most significantly differentially expressed genes reported in another large case-control study of MDD (
      • Mostafavi S.
      • Battle A.
      • Zhu X.
      • Potash J.B.
      • Weissman M.M.
      • Shi J.
      • et al.
      Type I interferon signaling genes in recurrent major depression: Increased expression detected by whole-blood RNA sequencing.
      ); only SRSF5 was underexpressed in both of these lists.
      To assess the robustness of our results to key modeling assumptions, we conducted two sensitivity analyses. First, we repeated the analysis using a model for differential expression that did not include comorbid anxiety as an explanatory factor. Only 24 genes were identified by Bayesian analysis as significantly coexpressed; of these 24 genes, 21 were concordant for sign of fold change, and 15 of these 21 genes were also included in the MDD-165 consensus set (Supplement). The reduced sensitivity when not adjusting for anxiety likely reflects the importance of refining or homogenizing the clinical phenotype in the search for biomarkers. Second, we repeated the analysis using a model for differential expression that included BMI as a covariate. BMI was not measured in the GSK–HiTDiP study, so this analysis was restricted to 74 patients and 80 control subjects from the Janssen–BRC dataset. Inclusion of BMI in the model tended to increase the p values for MDD-related differences, and only 73 genes in the MDD-165 consensus set were differentially expressed at the probability threshold q2 = .0246 (Supplemental Figure S3). The reduced sensitivity of the analysis including BMI as a covariate likely reflects the known proinflammatory effect of obesity and the reduced sample size available.

      Gene Set Enrichment Analyses Across GO Terms

      We conducted a series of enrichment analyses to functionally characterize the MDD-165 consensus set. The top 10 most significantly enriched GO terms all were related to the immune response to infection (Bonferroni q < .05 for the top 9 terms) (Supplemental Table S7). Among the MDD-165 consensus set, 90 genes were overexpressed in cases compared with controls. The UP gene list was significantly enriched for GO terms related to the innate immune system and response to infection (Bonferroni q < .05 for the top 17 terms) (Table 1). In addition, 75 genes were significantly underexpressed in cases compared with controls. As we hypothesized based on prior literature (
      • Blume J.
      • Douglas S.D.
      • Evans D.L.
      Immune suppression and immune activation in depression.
      ), this gene list was enriched for GO terms related to the adaptive immune system (p ≤ 10−3 for top 20 GO terms; Table 1), but there was no significant enrichment of the DOWN genes at the more stringent Bonferroni threshold.
      Table 1Enrichment Analysis of the MDD-165 Consensus Set of Genes That Were Differentially Overexpressed (UP) or Underexpressed (DOWN) in Patients With MDD Compared With Healthy Control Subjects in Both GSK–HiTDiP and Janssen–BRC Studies
      RankGO:IDTermAnnotatedFoundExpectedFisher’s Test
      Overexpressed UP (Innate Immunity)
      1GO:0009617Response to bacterium405192.297.7e-13
      2GO:0042742Defense response to bacterium164130.937.0e-12
      3GO:0043207Response to external biotic stimulus651203.693.7e-10
      4GO:0051707Response to other organism651203.693.7e-10
      5GO:0009607Response to biotic stimulus680203.858.0e-10
      6GO:0002376Immune system process21493512.179.8e-10
      7GO:0006955Immune response1310267.426.3e-09
      8GO:0098542Defense response to other organism351141.999.1e-09
      9GO:0051704Multi-organism process21763212.329.2e-08
      10GO:0050832Defense response to fungus2150.129.7e-08
      11GO:0032496Response to lipopolysaccharide253111.431.7e-07
      12GO:0002237Response to molecule of bacterial origin265111.52.8e-07
      13GO:0006952Defense response1372247.773.2e-07
      14GO:0009605Response to external stimulus19472911.023.8e-07
      15GO:0045087Innate immune response820184.645.5e-07
      16GO:0009620Response to fungus3450.191.3e-06
      17GO:0019731Antibacterial humoral response2340.137.8e-06
      18GO:0031640Killing of cells of other organism2540.141.1e-05
      19GO:0044364Disruption of cells of other organism2540.141.1e-05
      20GO:0019730Antimicrobial humoral response2740.151.5e-05
      Underexpressed DOWN (Adaptive Immunity)
      1GO:0050851Antigen receptor-mediated signaling path13350.590.00029
      2GO:0046649Lymphocyte activation53292.340.0005
      3GO:0009059Macromolecule biosynthetic process43363219.090.00052
      4GO:0034645Cellular macromolecule biosynthetic proc41833118.420.00064
      5GO:0010557Positive regulation of macromolecule bio1361155.990.00066
      6GO:0032774RNA biosynthetic process32752614.420.00086
      7GO:0046631Alpha-beta T cell activation9940.440.00094
      8GO:0016070RNA metabolic process38732917.060.00095
      9GO:0050776Regulation of immune response721103.180.00112
      10GO:0006351Transcription, DNA-templated31612513.920.00122
      11GO:0097659Nucleic acid-templated transcription31772513.990.00131
      12GO:0042110T cell activation38271.680.0014
      13GO:0070489T cell aggregation38271.680.0014
      14GO:0071593Lymphocyte aggregation38471.690.00144
      15GO:0031328Positive regulation of cellular biosyn1470156.470.00147
      16GO:0002562Somatic diversification of immune recept5230.230.00154
      17GO:0016444Somatic cell DNA recombination5230.230.00154
      18GO:0070486Leukocyte aggregation39071.720.00158
      19GO:0045893Positive regulation of transcription, DN1178135.190.00161
      20GO:1903508Positive regulation of nucleic acid-temp1178135.190.00161
      The top 20 GO terms for biological processes are ranked according to their p values by Fisher’s exact test for significant enrichment. Bonferroni correction specifies a p-value threshold of p = 4.94 × 10−6 to achieve significance at q < 0.05. Such a correction is well known to be too stringent given that many of the GO terms are correlated; however, the top 17 terms were significantly enriched in the overexpressed UP gene set after Bonferroni correction.
      GO, Gene Ontology; GSK–HiTDiP, GlaxoSmithKline–High-Throughput Disease-specific target Identification Program; Janssen–BRC, Janssen–Brain Resource Company; MDD, major depressive disorder; MDD-165 consensus set, set of 165 replicably concordant genes.

      Protein Interaction and Whole-Genome Transcriptional Network Analyses

      We used Search Tool for the Retrieval of Interacting Genes/Proteins analysis (
      • Szklarczyk D.
      • Franceschini A.
      • Wyder S.
      • Forslund K.
      • Heller D.
      • Huerta-Cepas J.
      • et al.
      STRING v10: Protein–protein interaction networks, integrated over the tree of life.
      ) to visualize the network of protein–protein interactions between the proteins coded by the MDD-165 consensus set. These genes were significantly enriched for protein–protein interactions (permutation test, p = 7 × 10−4) that were concentrated around mitogen-activated protein kinase 14 (MAPK14) and matrix metalloproteinase 9 (MMP9), which thus can be regarded as the most highly interactive hubs of this immune signaling network (Figure 1). Many of the proteins in this network are coded by genes that were independently identified by the NESDA study as differentially expressed in MDD; see Figure 1.
      Figure thumbnail gr1
      Figure 1Protein–protein interaction network for proteins coded by the set of 165 replicably concordant genes differentially expressed in both case-control studies of major depressive disorder. The network is represented by an undirected graph where links correspond to known protein–protein interactions and weights are proportional to the Search Tool for the Retrieval of Interacting Genes/Proteins confidence score
      (
      • Franceschini A.
      • Szklarczyk D.
      • Frankild S.
      • Kuhn M.
      • Simonovic M.
      • Roth A.
      • et al.
      STRING v9: 1. Protein–protein interaction networks, with increased coverage and integration.
      )
      . Only high-confidence (>0.7) links are retained, and disconnected genes are not shown. Red (blue) nodes correspond to genes over- (under-) expressed in major depressive disorder in both the GlaxoSmithKline–High-Throughput Disease-specific target Identification Program and Janssen–Brain Resource Company datasets. Smaller inner circles highlight proteins that are coded by genes also differentially expressed and with the same sign of fold change in a third large independent case-control study of major depressive disorder (Netherlands Study of Depression and Anxiety [NESDA])
      (
      • Jansen R.
      • Penninx B.W.J.H.
      • Madar V.
      • Xia K.
      • Milaneschi Y.
      • Hottenga J.J.
      • et al.
      Gene expression in major depressive disorder.
      )
      thresholded to control false discovery rate (FDR) at 20% (white circles), 10% (gray circles), and 5% (black circles).
      We constructed a graph representing significant coexpression of a pair of genes (nodes) as an edge drawn between them. As previously reported (
      • Langfelder P.
      • Horvath S.
      WGCNA: An R package for weighted correlation network analysis.
      ), this whole-genome transcriptional network or transcriptome had complex topological properties, including a community structure comprising modules of coexpressed genes enriched for specific GO terms (Figure 2). The UP genes were concentrated in the normative module (red) significantly enriched for innate immune response GO terms (e.g., myeloid cell activation involved in immune response, q = 5 × 10−5, Bonferroni corrected), whereas the DOWN genes were concentrated in a normative module (pink) specialized for translation-based terms (Figure 2). We confirmed in a hypothesis-driven analysis that the module associated with DOWN gene expression was also significantly enriched for adaptive immune response terms (e.g., T-cell differentiation, p = .0002). There were significant case-control differences in eigengene scores summarizing expression of all genes within modules enriched for innate immune response (red, p = .031; yellow, p = .015), translation (pink, p = .002), and one additional module with no significant enrichment terms (tan, p = .025); see Supplemental Table S11.
      Figure thumbnail gr2
      Figure 2Major depressive disorder (MDD)-related genes in the context of the normative whole-genome transcriptome. Overexpressed genes (or UP genes) in patients with MDD are concentrated in a module of the normative gene coexpression network specialized for innate immune response, whereas underexpressed genes (or DOWN genes) are concentrated in a module partially specialized for adaptive immune response. (A) The modules of the normative transcriptome are highlighted in different colors. (B) The MDD-related genes are colored according to their normative module affiliation, and representative genes are text-labeled. (C) The MDD-related genes are colored green for overexpressed (or UP) genes and are colored red for underexpressed (or DOWN) genes. The text labels highlight the functions of the corresponding modules of the normative transcriptome. rRNA, ribosomal RNA; SRP, signal recognition particle.

      Cell Type–Specific Expression Patterns for the MDD-165 Consensus Set

      To assess the cellular specificity of the MDD-165 consensus set, we used microarray data from an independent study of cell type–specific gene expression across eight major classes of immune cells in healthy control subjects (
      • Mabbott N.A.
      • Baillie J.K.
      • Brown H.
      • Freeman T.C.
      • Hume D.A.
      An expression atlas of human primary cells: Inference of gene function from coexpression networks.
      ). Remarkably, 89 of the 90 UP genes formed four clusters of gene coexpression, three of which represented strong overexpression of a subset of UP genes in one or two myeloid cell classes (Figure 3). Likewise, 71 of the 75 DOWN genes formed two clusters of gene coexpression, each representing cell-specific overexpression of a subset of DOWN genes in one or two lymphoid or erythroblast classes (Figure 3).
      Figure thumbnail gr3
      Figure 3Cell class–specific expression patterns for the set of 165 replicably concordant genes. We estimated the correlation between each possible pair of the 165 major depressive disorder (MDD)-related genes and identified clusters of genes with similar expression patterns in an independent microarray dataset on specific cell types. The set of 165 replicably concordant genes formed six clusters, with each cluster comprising a subset of genes that had strong mutual coexpression across a range of eight distinct cell classes: erythroblasts, monocytes, monocyte-derived dendritic cells (MDDCs), neutrophils, B cells, CD4+ T cells, CD8+ T cells, and natural killer (NK) cells. (A) The six clusters of genes with strongly correlated expression profiles. Clusters 1 to 4 (left) comprised genes that were overexpressed in MDD (or UP genes), and clusters 5 and 6 (right) comprised genes that were underexpressed in MDD (or DOWN genes). (B) Histograms of clustered gene expression across cell types for each of the four clusters of UP genes overexpressed in MDD (from top to bottom: clusters 1–4). The x-axis color legend codes for different cell types. (C) Histograms of clustered gene expression across cell types for each of the two clusters of DOWN genes underexpressed in MDD (from top to bottom: clusters 5 and 6). The x-axis color legend codes for different cell types.

      Coupled UP and DOWN Gene Expression in Patients With MDD

      To explore the extent to which UP gene overexpression was related to DOWN gene underexpression, we estimated the mean UP gene expression (over all 90 differentially overexpressed genes) and the mean DOWN gene expression (over all 75 differentially underexpressed genes) for each of 113 patients with MDD in the GSK–HiTDiP study and each of 94 patients with MDD in the Janssen–BRC study. In both studies, UP and DOWN expression scores were strongly negatively correlated (Figure 4), indicating that these case-control differences were coupled at the level of individual patients (Janssen–BRC: r = −.82, p < 10−4; GSK–HiTDiP: r = −.74, p < 10−4).
      Figure thumbnail gr4
      Figure 4Opposing and coordinated expression of innate and adaptive immune transcripts in patients with major depressive disorder (MDD). (A, B) Scatter plots for the mean expression across the 75 underexpressed (DOWN) genes plotted against the mean expression of 90 overexpressed (UP) genes in the Janssen–Brain Resource Company (Janssen–BRC) (A) and GlaxoSmithKline–High-Throughput Disease-specific target Identification Program (GSK–HiTDiP) (B) datasets. Each point corresponds to a patient. Blue indicates case, and gray indicates control (CTL). Regression lines are shown in blue and gray, respectively. In (A), an individual patient’s data point (red outline) is projected onto the regression line (red circle). The distance from the origin to the point on the regression line is the bioscalar value for that patient. Inset illustrates the projection (black) of all individual patient data points (blue) onto the sample regression line. (C, D) Left panels: Box plot of MDD-165 bioscalar values in controls and cases. Green line indicates the threshold identifying the top third (tertile) of the MDD-165 bioscalar distribution as a subgroup of inflamed patients with MDD (designated immune-MDD). Right panels: Receiver operating characteristic curves were calculated for MDD vs. CTL and immune-MDD vs. CTL classifications. (E) Correlation of each depressed patient’s MDD-165 bioscalar with body mass index (BMI). (F, G) Box plots of MDD-165 bioscalar for subgroups defined by a comorbid diagnosis of anxiety disorder. AUC, area under the curve.
      Given the strong linear relationship between UP and DOWN gene expression in patients with MDD, we mapped each patient’s data to a point on the fitted regression line, thereby characterizing each patient by a bioscalar that locates it on a single axis or dimension of coupled UP/DOWN gene expression; see Figure 4A (inset). As expected, there were significant case-control differences in the MDD-165 bioscalar in both studies (Cohen’s d = 0.68 for GSK–HiTDiP and 0.58 for Janssen–BRC, corresponding to a medium effect size; Figure 4C, D), with patients on average having more positive values, indicating a greater shift in the direction of coupled innate activation and adaptive inactivation. Receiver operating characteristic analysis indicated that cases and controls were classified with an area under the curve of 0.71 for the GSK–HiTDiP study and an area under the curve of 0.67 for the Janssen–BRC study (Figure 4C, D). This moderately accurate classification performance reflects the fact that the bioscalar distributions of cases and controls are overlapping. We also assessed the performance of the bioscalar not to discriminate cases from controls (which is ultimately a clinical diagnostic decision) but rather to identify the top third “most inflamed” patients with MDD. To do this, we defined a cutoff value for the bioscalar corresponding to the top tertile critical value of its distribution in each study (0.57). Participants with MDD-165 > 0.57 were classified correctly as belonging to such an inflamed MDD subgroup (sensitivity = 100% by definition) with a high specificity (93% and 88% in GSK–HiTDiP and Janssen–BRC studies, respectively).
      We provisionally explored correlations between the MDD-165 bioscalar and sociodemographic or clinical differences among the patients with MDD in the Janssen–BRC study. BMI (r72 = .20, p = .099), the CORE total score for melancholic symptom severity (r88 = .18, p = .097), and the presence or absence of comorbid substance abuse disorder (two-tailed t63 = 1.9, p = .07) were not significantly associated with bioscalar scores. There was a significant difference in bioscalar scores between subgroups of patients with MDD with or without comorbid anxiety disorder (t83 = −2.4, p = .02).

      Discussion

      We have reported the differential expression analysis of whole-genome microarray data from two independent case-control studies of patients with MDD compared with healthy control subjects. Using Bayesian methods, we identified a set of 165 consensus genes that were replicably associated with MDD, sharing the same direction of fold change in the cases compared with controls in both studies. The robustness of these results was further supported by comparison with the results of the largest prior case-control microarray study of depression, the NESDA study (
      • Jansen R.
      • Penninx B.W.J.H.
      • Madar V.
      • Xia K.
      • Milaneschi Y.
      • Hottenga J.J.
      • et al.
      Gene expression in major depressive disorder.
      ); several genes differentially expressed in the NESDA study were likewise differentially expressed in these data, for example, MAPK14 and MMP9. Both of these genes were consistently overexpressed in patients with MDD and both code proteins (MAPK14 and matrix metallopeptidase 9) that were hubs of a network of interactions between immune signaling proteins coded by many of the other differentially expressed genes. A MAPK14 inhibitor has been tried for treatment of major depression but did not demonstrate consistently significant effects on symptom rating scales compared with placebo at the single dose tested (
      • Inamdar A.
      • Merlo-Pich E.
      • Gee M.
      • Makumi C.
      • Mistry P.
      • Robertson J.
      • et al.
      Evaluation of antidepressant properties of the p38 MAP kinase inhibitor losmapimod (GW856553) in major depressive disorder: Results from two randomised, placebo-controlled, double-blind, multicentre studies using a Bayesian approach.
      ).
      The consensus set of 165 genes associated with MDD was divided into two approximately equal-sized subsets; here, 90 so-called UP genes were overexpressed in patients and 75 so-called DOWN genes were underexpressed in patients. The overexpressed UP genes were significantly enriched for GO terms related to the response to infection and the innate immune system. These gene transcripts were not functionally unrelated to or independent of each other. Most UP genes were affiliated with the module of the normative gene transcriptional network specialized for innate immune response. Smaller clusters of highly correlated UP genes were typically enriched in neutrophils, monocytes, and monocyte-derived dendritic cells. These results are compatible with prior data implicating activation of the innate immune system and increased proinflammatory cytokine signaling in (some people with) depression (
      • Dantzer R.
      • O’Connor J.
      • Freund G.
      • Johnson R.
      • Kelley K.
      From inflammation to sickness and depression: When the immune system subjugates the brain.
      ,
      • Miller A.
      • Maletic V.
      • Raison C.
      Inflammation and its discontents: The role of cytokines in the pathophysiology of major depression.
      ).
      In contrast, we found that the underexpressed DOWN genes were enriched for GO terms related to T-cell function and adaptive immunity. The DOWN genes were affiliated with modules of the transcriptome partially specialized for adaptive immune function, and smaller clusters of strongly coexpressed DOWN genes were enriched in T cells, B cells, and NK cells. These results are compatible with prior data implicating relative suppression of the cellular immune system in (some people with) depression, including evidence for decreased NK cell cytotoxicity and decreased proliferation of lymphocytes challenged with mitogens in vitro (
      • Zorrilla E.P.
      • Luborsky L.
      • McKay J.R.
      • Rosenthal R.
      • Houldin A.
      • Tax A.
      • et al.
      The relationship of depression and stressors to immunological assays: A meta-analytic review.
      ), and with recent data suggesting that deficiencies of NK and T cells co-occur with inflammatory monocyte activation as related phenomena in the same patients with MDD (
      • Grosse L.
      • Hoogenboezem T.
      • Ambrée O.
      • Bellingrath S.
      • Jörgens S.
      • de Wit H.J.
      • et al.
      Deficiencies of the T and natural killer cell system in major depressive disorder: T regulatory cell defects are associated with inflammatory monocyte activation.
      ,
      • Snijders G.
      • Schiweck C.
      • Mesman E.
      • Grosse L.
      • de Wit H.J.
      • Nolen W.A.
      • et al.
      A dynamic course of T cell defects in individuals at risk for mood disorders.
      ,
      • Carvalho L.A.
      • Bergink V.
      • Sumaski L.
      • Wijkhuijs J.
      • Hoogendijk W.J.
      • Birkenhager T.K.
      • et al.
      Inflammatory activation is associated with a reduced glucocorticoid receptor alpha/beta expression ratio in monocytes of inpatients with melancholic major depressive disorder.
      ,
      • Weigelt K.
      • Carvalho L.A.
      • Drexhage R.C.
      • Wijkhuijs A.
      • de Wit H.
      • van Beveren N.J.
      • et al.
      TREM-1 and DAP12 expression in monocytes of patients with severe psychiatric disorders: EGR3, ATF3 and PU.1 as important transcription factors.
      ).
      In both studies, we also observed a strong statistical association between mean UP gene overexpression and mean DOWN gene underexpression in each patient. This provided robust evidence that these complementary immunophenotypes are indeed coupled at the level of an individual patient. It also motivated the concept that each individual patient may be located by a single number (bioscalar) on a spectrum of coupled change in innate and adaptive immune system function. Patients with MDD had significantly higher scores on this bioscalar, indicating a mean shift to relatively increased UP gene expression and decreased DOWN gene expression. The subgroup of most-inflamed patients with MDD, defined as the top tertile of the MDD-165 bioscalar distribution, were identified with high sensitivity and specificity, suggesting that this may in the future prove to be a useful biomarker for defining an abnormally inflamed subgroup of patients with MDD.
      So far, we have used the terms overexpression and underexpression simply to describe higher and lower levels of measured mRNA in cases compared with controls. At least three (mutually nonexclusive) explanatory factors are plausible: cellular, genetic, and environmental. First, case-control differences in expression could reflect differences in cell counts. For example, increased monocyte counts have been reported in MDD (
      • Maes M.
      • Lambrechts J.
      • Bosmans E.
      • Jacobs J.
      • Suy E.
      • Vandervorst C.
      • et al.
      Evidence for a systemic immune activation during depression: Results of leukocyte enumeration by flow cytometry in conjunction with monoclonal antibody staining.
      ) and could cause apparent overexpression of innate immune system genes measured in whole-blood samples (and relative underexpression of adaptive immune genes). It will be important in the future to combine cytometry and transcriptional measurements in the same patients and to measure case-control expression differences in sorted cell samples. Second, mRNA changes could be quantitative traits determined by DNA variation at expression quantitative trait loci. This genetic explanation would require that (at least some) patients with MDD had a consistent profile of allelic variation in inflammation-related genes (
      • Wong M.L.
      • Dong C.
      • Maestre-Mesa J.
      • Licinio J.
      Polymorphisms in inflammation-related genes are associated with susceptibility to major depression and antidepressant response.
      ), but large genome-wide association studies have so far failed to identify genetic variants robustly linked to risk for major depression (
      • Ripke S.
      • Wray N.R.
      • Lewis C.M.
      • Hamilton S.P.
      • Weissman M.M.
      • et al.
      Major Depressive Disorder Working Group of the Psychiatric GWAS Consortium
      A mega-analysis of genome-wide association studies for major depressive disorder.
      ). The absence of significant genome-wide association study findings in MDD could be regarded as problematic for an expression quantitative trait locus interpretation of expression changes, or it could be discounted on the grounds that classically designed association studies were underpowered to detect DNA variations occurring in only a subgroup of immunologically dysfunctional patients. Third, gene expression changes could have been induced in patients by exposure to some shared environmental stimulus. The risk for depression is associated with adverse events in the biological and social environments such as infection, childhood abuse, and bereavement (
      • Danese A.
      • Pariante C.M.
      • Caspi A.
      • Taylor A.
      • Poulton R.
      Childhood maltreatment predicts adult inflammation in a life-course study.
      ). In addition, psychosocial stress has been linked to peripheral immune state changes, such as increased proinflammatory cytokines and monocyte activation (
      • Rohleder N.
      Stimulation of systemic low-grade inflammation by psychosocial stress.
      ), that are compatible with overexpression of innate immune system genes.
      It is a methodological limitation that the GSK–HiTDiP study, designed primarily to identify risk genes for MDD, did not provide data on severity of depressive symptoms or BMI. The case-control comparisons are not controlled for potentially confounding effects of cigarette smoking, race, comorbid medical disease, or socioeconomic status on peripheral immune status. Future studies of transcriptional biomarkers in MDD, with more detailed clinical and immunological phenotyping and more complete control of potential confounds, will be required to evaluate the generalizability of these results.
      In short, we have reported replicable new evidence in support of peripheral immune gene expression markers for MDD. It may be fruitful in the future to investigate coupled overexpression of innate immune genes and underexpression of adaptive immune genes as a predictor of antidepressant response to novel immunotherapeutics.

      Acknowledgments and Disclosures

      This work was supported by a grant from the Medical Research Council to the Immunopsychiatry Consortium led by ETB (Grant No. MR/L014815/1). Data collection for the GSK–HiTDiP case-control study was funded by GlaxoSmithKline; data collection for the Janssen–BRC study was funded by Janssen. PEV is supported by a Medical Research Council Bioinformatics Research Fellowship (Grant No. MR/K020706/1). TCF is supported by an Institute Strategic Grant from the Biotechnology and Biological Sciences Research Council (Grant No. BB/JO1446X/1). CMP is supported by the National Institute for Health Research (NIHR) Mental Health Biomedical Research Centre at South London and Maudsley NHS Foundation Trust and King’s College London. ETB is supported by the NIHR Cambridge Biomedical Research Centre.
      We thank David Hume for comments on the manuscript; Fiona Kelly and Stewart Bates for supporting access to the GSK–HiTDiP dataset; Pieter Peeters, Wilhelmus Drinkenberg, and William Talloan for supporting access to the Janssen–BRC and HiTDiP datasets; and Evian Gordon from the BRC.
      The microarray and clinical data from the Janssen-BRC study are accessible to scientists through the independent YODA Project process (http://yoda.yale.edu). The GSK-HiTDiP data are accessible on GEO (https://www.ncbi.nlm.nih.gov/geo) as GSE98793. For more information on combining the datasets please contact the YODA Project ( [email protected] ).
      GMW and WCD are employees and stockholders of Janssen Research & Development of Johnson & Johnson. SK and RH are employees and stockholders of GlaxoSmithKline (GSK). ETB is employed half-time by GSK and half-time by the University of Cambridge; he holds stock in GSK. ETB is also a shareholder in the BRC. CMP has received research funding from Johnson & Johnson and speaker fees from Lundbeck. The remaining authors all report no biomedical financial interests or potential conflicts of interest.

      Supplementary Material

      References

        • Wium-Andersen M.K.
        • Ørsted D.D.
        • Nielsen S.F.
        • Nordestgaard B.G.
        Elevated C-reactive protein levels, psychological distress, and depression in 73,131 individuals.
        JAMA Psychiatry. 2013; 70: 176-184
        • Valkanova V.
        • Ebmeier K.P.
        • Allan C.L.
        CRP, IL-6 and depression: A systematic review and meta-analysis of longitudinal studies.
        J Affect Disord. 2013; 150: 736-744
        • Haapakoski R.
        • Mathieu J.
        • Ebmeier K.P.
        • Alenius H.
        • Kivimaki M.
        Cumulative meta-analysis of interleukins 6 and 1β, tumour necrosis factor α and C-reactive protein in patients with major depressive disorder.
        Brain Behav Immun. 2015; 49: 206-215
        • Matcham F.
        • Rayner L.
        • Steer S.
        • Hotopf M.
        The prevalence of depression in rheumatoid arthritis: A systematic review and meta-analysis.
        Rheumatology. 2013; 52: 2136-2148
        • Khandaker G.M.
        • Pearson R.M.
        • Zammit S.
        • Lewis G.
        • Jones P.B.
        Association of serum interleukin 6 and C-reactive protein in childhood with depression and psychosis in young adult life: A population-based longitudinal study.
        JAMA Psychiatry. 2014; 71: 1121-1128
        • Bull S.J.
        • Huezo-Diaz P.
        • Binder E.B.
        • Cubells J.F.
        • Ranjith G.
        • Maddock C.
        • et al.
        Functional polymorphisms in the interleukin-6 and serotonin transporter genes, and depression and fatigue induced by interferon-α and ribavirin treatment.
        Mol Psychiatry. 2009; 14: 1095-1104
        • Hepgul N.
        • Cattaneo A.
        • Agarwal K.
        • Baraldi S.
        • Borsini A.
        • Bufalino C.
        • et al.
        Transcriptomics in interferon-α-treated patients identifies inflammation-, neuroplasticity- and oxidative stress-related signatures as predictors and correlates of depression.
        Neuropsychopharmacology. 2016; 41: 2502-2511
        • Dantzer R.
        • O’Connor J.
        • Freund G.
        • Johnson R.
        • Kelley K.
        From inflammation to sickness and depression: When the immune system subjugates the brain.
        Nat Rev Neurosci. 2008; 9: 46-56
        • Banks W.A.
        • Kastin A.J.
        • Gutierrez E.G.
        Penetration of interleukin-6 across the murine blood-brain barrier.
        Neurosci Lett. 1994; 179: 53-56
        • Perry V.H.
        • Cunningham C.
        • Holmes C.
        Systemic infections and inflammation affect chronic neurodegeneration.
        Nat Rev Immunol. 2007; 7: 161-167
        • Miller A.H.
        • Raison C.L.
        Are anti-inflammatory therapies viable treatments for psychiatric disorders? Where the rubber meets the road.
        JAMA Psychiatry. 2015; 72: 527-528
        • Raison C.
        • Capuron L.
        • Miller A.
        Cytokines sing the blues: Inflammation and the pathogenesis of depression.
        Trends Immunol. 2006; 27: 24-31
        • Lin E.
        • Tsai S.J.
        Genome-wide microarray analysis of gene expression profiling in major depression and antidepressant therapy.
        Prog Neuropsychopharmacol Biol Psychiatry. 2016; 64: 334-340
        • Jansen R.
        • Penninx B.W.J.H.
        • Madar V.
        • Xia K.
        • Milaneschi Y.
        • Hottenga J.J.
        • et al.
        Gene expression in major depressive disorder.
        Mol Psychiatry. 2016; 21: 444
        • Garbett K.A.
        • Vereczkei A.
        • Kálmán S.
        • Brown J.A.
        • Taylor W.D.
        • Faludi G.
        • et al.
        Coordinated messenger RNA/microRNA changes in fibroblasts of patients with major depression.
        Biol Psychiatry. 2015; 77: 256-265
        • Mostafavi S.
        • Battle A.
        • Zhu X.
        • Potash J.B.
        • Weissman M.M.
        • Shi J.
        • et al.
        Type I interferon signaling genes in recurrent major depression: Increased expression detected by whole-blood RNA sequencing.
        Mol Psychiatry. 2014; 19: 1267-1274
        • Glahn D.C.
        • Curran J.E.
        • Winkler A.M.
        • Carless M.A.
        • Kent J.W.
        • Charlesworth J.C.
        • et al.
        High dimensional endophenotype ranking in the search for major depression risk genes.
        Biol Psychiatry. 2012; 71: 6-14
        • Shelton R.C.
        • Claiborne J.
        • Sidoryk-Wegrzynowicz M.
        • Reddy R.
        • Aschner M.
        • Lewis D.A.
        • et al.
        Altered expression of genes involved in inflammation and apoptosis in frontal cortex in major depression.
        Mol Psychiatry. 2011; 16: 751-762
        • Blume J.
        • Douglas S.D.
        • Evans D.L.
        Immune suppression and immune activation in depression.
        Brain Behav Immun. 2011; 25: 221-229
        • Miller A.
        • Maletic V.
        • Raison C.
        Inflammation and its discontents: The role of cytokines in the pathophysiology of major depression.
        Biol Psychiatry. 2009; 65: 732-741
        • Grosse L.
        • Hoogenboezem T.
        • Ambrée O.
        • Bellingrath S.
        • Jörgens S.
        • de Wit H.J.
        • et al.
        Deficiencies of the T and natural killer cell system in major depressive disorder: T regulatory cell defects are associated with inflammatory monocyte activation.
        Brain Behav Immun. 2016; 54: 38-44
        • Snijders G.
        • Schiweck C.
        • Mesman E.
        • Grosse L.
        • de Wit H.J.
        • Nolen W.A.
        • et al.
        A dynamic course of T cell defects in individuals at risk for mood disorders.
        Brain Behav Immun. 2016; 58: 11-17
        • Grosse L.
        • Carvalho L.A.
        • Wijkhuijs A.M.
        • Bellingrath S.
        • Ruland T.
        • Ambree O.
        • et al.
        Clinical characteristics of inflammation-associated depression: Monocyte gene expression is age-related in major depressive disorder.
        Brain Behav Immun. 2015; 44: 48-56
        • Domenici E.
        • Willé D.R.
        • Tozzi F.
        • Prokopenko I.
        • Miller S.
        • McKeown A.
        • et al.
        Plasma protein biomarkers for depression and schizophrenia by multi analyte profiling of case-control collections.
        PLoS One. 2010; 5: e9166
        • Muglia P.
        • Tozzi F.
        • Galwey N.W.
        • Francks C.
        • Upmanyu R.
        • Kong X.Q.
        • et al.
        Genome-wide association study of recurrent major depressive disorder in two European case-control cohorts.
        Mol Psychiatry. 2010; 15: 589-601
        • Tilahun A.
        • Lin D.
        • Shkedy Z.
        • Geys H.
        • Alonso A.
        • Peeters P.
        • et al.
        Genomic biomarkers for depression: Feature-specific and joint biomarkers.
        Stat Biopharm Res. 2010; 2: 419-434
        • Zoon H.
        • Veth C.
        • Arns M.
        • Drinkenburg W.
        • Talloen W.
        • Peeters P.
        • et al.
        EEG alpha power as an intermediate measure between brain-derived neurotrophic factor Val66Met and depression severity in patients with major depressive disorder.
        J Clin Neurophys. 2013; 30: 261-267
        • Liu Y.
        • Yieh L.
        • Yang T.
        • Drinkenburg W.
        • Peeters P.
        • Steckler T.
        • et al.
        Metabolomic biosignature differentiates melancholic depressive patients from healthy controls.
        BMC Genomics. 2016; 17: 669
        • Sheehan D.V.
        • Lecrubier Y.
        • Sheehan K.H.
        • Amorim P.
        • Janavs J.
        • Weiller E.
        • et al.
        The Mini-International Neuropsychiatric Interview (MINI): The development and validation of a structured diagnostic psychiatric interview for DSM-IV and ICD-10.
        J Clin Psychiatry. 1998; 59: 22-33
        • Ritchie M.E.
        • Phipson B.
        • Wu D.
        • Hu Y.
        • Law C.W.
        • Shi W.
        • et al.
        limma powers differential expression analyses for RNA-sequencing and microarray studies.
        Nucleic Acids Res. 2015; 43: e47
        • Smyth G.K.
        Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.
        Stat Appl Genet Mol Biol. 2004; 3: 3
        • Blangiardo M.
        • Richardson S.
        Statistical tools for synthesizing lists of differentially expressed features in related experiments.
        Genome Biol. 2007; 8: R54
        • Blangiardo M.
        • Cassese A.
        • Richardson S.
        sdef: An R package to synthesize lists of significant features in related experiments.
        BMC Bioinformatics. 2010; 11: 270
        • Zhang B.
        • Horvath S.
        A general framework for weighted gene co-expression network analysis.
        Stat Appl Genet Mol Biol. 2005; 4: 17
      1. Alexa A, Rahnenfuhrer J (2010): topGO: Enrichment analysis for gene ontology. R package, version 2.0.

        • Ashburner M.
        • Ball C.A.
        • Blake J.A.
        • Botstein D.
        • Butler H.
        • Cherry J.M.
        • et al.
        Gene Ontology: Tool for the unification of biology.
        Nat Genet. 2000; 25: 25-29
        • Szklarczyk D.
        • Franceschini A.
        • Wyder S.
        • Forslund K.
        • Heller D.
        • Huerta-Cepas J.
        • et al.
        STRING v10: Protein–protein interaction networks, integrated over the tree of life.
        Nucleic Acids Res. 2015; 43: D447-D452
        • Franceschini A.
        • Szklarczyk D.
        • Frankild S.
        • Kuhn M.
        • Simonovic M.
        • Roth A.
        • et al.
        STRING v9: 1. Protein–protein interaction networks, with increased coverage and integration.
        Nucleic Acids Res. 2013; 41: D808-D815
        • Mabbott N.A.
        • Baillie J.K.
        • Brown H.
        • Freeman T.C.
        • Hume D.A.
        An expression atlas of human primary cells: Inference of gene function from coexpression networks.
        BMC Genomics. 2013; 14: 632
        • Freeman T.C.
        • Goldovsky L.
        • Brosch M.
        • Van Dongen S.
        • Mazière P.
        • Grocock R.J.
        • et al.
        Construction, visualisation, and clustering of transcription networks from microarray expression data.
        PLoS Comput Biol. 2007; 3: e206
        • Theocharidis A.
        • van Dongen S.
        • Enright A.
        • Freeman T.
        Network visualization and analysis of gene expression data using BioLayout Express3D.
        Nat Protoc. 2009; 4: 1535-1550
        • Langfelder P.
        • Horvath S.
        WGCNA: An R package for weighted correlation network analysis.
        BMC Bioinformatics. 2008; 9: 559
        • Inamdar A.
        • Merlo-Pich E.
        • Gee M.
        • Makumi C.
        • Mistry P.
        • Robertson J.
        • et al.
        Evaluation of antidepressant properties of the p38 MAP kinase inhibitor losmapimod (GW856553) in major depressive disorder: Results from two randomised, placebo-controlled, double-blind, multicentre studies using a Bayesian approach.
        J Psychopharmacol. 2014; 28: 570-581
        • Zorrilla E.P.
        • Luborsky L.
        • McKay J.R.
        • Rosenthal R.
        • Houldin A.
        • Tax A.
        • et al.
        The relationship of depression and stressors to immunological assays: A meta-analytic review.
        Brain Behav Immun. 2001; 15: 199-226
        • Carvalho L.A.
        • Bergink V.
        • Sumaski L.
        • Wijkhuijs J.
        • Hoogendijk W.J.
        • Birkenhager T.K.
        • et al.
        Inflammatory activation is associated with a reduced glucocorticoid receptor alpha/beta expression ratio in monocytes of inpatients with melancholic major depressive disorder.
        Transl Psychiatry. 2014; 4: e344
        • Weigelt K.
        • Carvalho L.A.
        • Drexhage R.C.
        • Wijkhuijs A.
        • de Wit H.
        • van Beveren N.J.
        • et al.
        TREM-1 and DAP12 expression in monocytes of patients with severe psychiatric disorders: EGR3, ATF3 and PU.1 as important transcription factors.
        Brain Behav Immun. 2011; 25: 1162-1169
        • Maes M.
        • Lambrechts J.
        • Bosmans E.
        • Jacobs J.
        • Suy E.
        • Vandervorst C.
        • et al.
        Evidence for a systemic immune activation during depression: Results of leukocyte enumeration by flow cytometry in conjunction with monoclonal antibody staining.
        Psychol Med. 1992; 22: 45-53
        • Wong M.L.
        • Dong C.
        • Maestre-Mesa J.
        • Licinio J.
        Polymorphisms in inflammation-related genes are associated with susceptibility to major depression and antidepressant response.
        Mol Psychiatry. 2008; 13: 800-812
        • Ripke S.
        • Wray N.R.
        • Lewis C.M.
        • Hamilton S.P.
        • Weissman M.M.
        • et al.
        • Major Depressive Disorder Working Group of the Psychiatric GWAS Consortium
        A mega-analysis of genome-wide association studies for major depressive disorder.
        Mol Psychiatry. 2013; 18: 497-511
        • Danese A.
        • Pariante C.M.
        • Caspi A.
        • Taylor A.
        • Poulton R.
        Childhood maltreatment predicts adult inflammation in a life-course study.
        Proc Natl Acad Sci U S A. 2007; 104: 1319-1324
        • Rohleder N.
        Stimulation of systemic low-grade inflammation by psychosocial stress.
        Psychosom Med. 2014; 76: 181-189