Advertisement

Delineating the Genetic Component of Gene Expression in Major Depression

  • Lorenza Dall’Aglio
    Affiliations
    Social, Genetic and Developmental Psychiatry Centre, Institute of Psychiatry, Psychology and Neuroscience, King's College London, London, United Kingdom

    Department of Child and Adolescent Psychiatry, Erasmus University Medical Center Rotterdam, Rotterdam, The Netherlands

    Generation R Study Group, Erasmus University Medical Center Rotterdam, Rotterdam, The Netherlands
    Search for articles by this author
  • Cathryn M. Lewis
    Affiliations
    Social, Genetic and Developmental Psychiatry Centre, Institute of Psychiatry, Psychology and Neuroscience, King's College London, London, United Kingdom

    Department of Medical and Molecular Genetics, Faculty of Life Sciences and Medicine, King’s College London, London, United Kingdom
    Search for articles by this author
  • Oliver Pain
    Correspondence
    Address correspondence to Oliver Pain, Ph.D.
    Affiliations
    Social, Genetic and Developmental Psychiatry Centre, Institute of Psychiatry, Psychology and Neuroscience, King's College London, London, United Kingdom
    Search for articles by this author
Open AccessPublished:September 12, 2020DOI:https://doi.org/10.1016/j.biopsych.2020.09.010

      Abstract

      Background

      Major depression (MD) is determined by a multitude of factors including genetic risk variants that regulate gene expression. We examined the genetic component of gene expression in MD by performing a transcriptome-wide association study (TWAS), inferring gene expression–trait relationships from genetic, transcriptomic, and phenotypic information.

      Methods

      Genes differentially expressed in depression were identified with the TWAS FUSION method, based on summary statistics from the largest genome-wide association analysis of MD (n = 135,458 cases, n = 344,901 controls) and gene expression levels from 21 tissue datasets (brain; blood; thyroid, adrenal, and pituitary glands). Follow-up analyses were performed to extensively characterize the identified associations: colocalization, conditional, and fine-mapping analyses together with TWAS-based pathway investigations.

      Results

      Transcriptome-wide significant differences between cases and controls were found at 94 genes, approximately half of which were novel. Of the 94 significant genes, 6 represented strong, colocalized, and potentially causal associations with depression. Such high-confidence associations include NEGR1, CTC-467M3.3, TMEM106B, LRFN5, ESR2, and PROX2. Lastly, TWAS-based enrichment analysis highlighted dysregulation of gene sets for, among others, neuronal and synaptic processes.

      Conclusions

      This study sheds further light on the genetic component of gene expression in depression by characterizing the identified associations, unraveling novel risk genes, and determining which associations are congruent with a causal model. These findings can be used as a resource for prioritizing and designing subsequent functional studies of MD.

      Keywords

      Major depression (MD) constitutes one of the largest contributors to global disability worldwide, affecting around 322 million people and accounting for approximately 50 million years lived with disability (
      World Health Organization
      Depression and Other Common Mental Disorders: Global Health Estimates..
      ). This disorder is of complex origin, being determined by the interplay of a multitude of environmental factors (e.g., life events) and genetic variations (
      American Psychiatric Association
      Diagnostic and Statistical Manual of Mental Disorders.
      ). The genetic contribution to MD has been shown by twin studies, which yielded heritability estimates of approximately 37% (
      • Sullivan P.F.
      • Neale M.C.
      • Kendler K.S.
      Genetic epidemiology of major depression: Review and meta-analysis.
      ), and genome-wide association studies (GWASs), which estimated a 9% heritability as captured by common single nucleotide polymorphisms (SNPs) (
      • Howard D.M.
      • Adams M.J.
      • Clarke T.K.
      • Hafferty J.D.
      • Gibson J.
      • Shirali M.
      • et al.
      Genome-wide meta-analysis of depression identifies 102 independent variants and highlights the importance of the prefrontal brain regions.
      ).
      GWASs have started to identify the specific genes underlying depression genetic risk. In one of the most recent and largest GWASs to date (n = 135,458 cases, n = 344,901 controls), Wray et al. (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ) uncovered 44 independent loci associated with MD. Drawing on a larger sample (n = 246,363 cases, n = 561,190 controls) and a broader phenotype definition, Howard et al. (
      • Howard D.M.
      • Adams M.J.
      • Clarke T.K.
      • Hafferty J.D.
      • Gibson J.
      • Shirali M.
      • et al.
      Genome-wide meta-analysis of depression identifies 102 independent variants and highlights the importance of the prefrontal brain regions.
      ) observed 102 independent loci, further demonstrating the polygenicity of the disorder. While GWASs have enabled the identification of SNPs conferring susceptibility to depression, the functional significance of such genetic variants remains to be elucidated. The examination of intermediary processes between genomics and the phenotype, such as gene expression, would permit greater insight into the molecular mechanisms underlying depression.
      Gene expression is a biological process permitting the translation of the genetic code into functional products, such as proteins and functional RNA. While gene expression is affected by extrinsic exposures, a substantial role for SNPs on the transcriptome has also been observed, with 50% of common variants being associated with gene expression, in any tissue (
      GTEx Consortium
      Genetic effects on gene expression across human tissues.
      ). Genetic variants affecting gene expression are referred to as expression quantitative trait loci (eQTL).
      Studies mapping eQTLs enabled the development of publicly available reference panels containing information on SNP-transcription relationships, in multiple tissues, across different samples. Based on these, gene expression can be predicted in any genetically mapped trait. The use of predicted gene expression allows for the analysis of transcription-trait associations with larger samples, inexpensively, and in multiple tissues, as opposed to classical observed gene expression studies. Importantly, such a novel approach solely captures the genetic component of gene expression, i.e., gene expression resulting from SNP variation (usually in cis, within 500 kb from gene boundaries). Differential expression that is a consequence of depression (reverse causation) can, therefore, be excluded.
      Several methods to predict transcription-trait associations are available. These leverage 3 key sources of information: 1) reference panels with data on SNPs-transcription associations (e.g., GTEx [Genotype-Tissue Expression] Consortium), 2) individual- or summary statistic–level GWAS data capturing SNP-trait associations, and 3) a linkage disequilibrium (LD) (i.e., correlations across SNPs) reference. Of the available methods, transcriptome-wide association studies (TWASs) (
      • Barbeira A.N.
      • Dickinson S.P.
      • Bonazzola R.
      • Zheng J.
      • Wheeler H.E.
      • Torres J.M.
      • et al.
      Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics.
      ,
      • Gusev A.
      • Ko A.
      • Shi H.
      • Bhatia G.
      • Chung W.
      • Penninx B.W.J.H.
      • et al.
      Integrative approaches for large-scale transcriptome-wide association studies.
      ) are the most commonly used. These explore the transcription-trait association with a gene-by-gene approach. The integration of SNP effects into genes lowers the multiple testing burden and increases power, thus allowing TWASs to identify genetic associations potentially overlooked by GWASs (
      • Gusev A.
      • Ko A.
      • Shi H.
      • Bhatia G.
      • Chung W.
      • Penninx B.W.J.H.
      • et al.
      Integrative approaches for large-scale transcriptome-wide association studies.
      ).
      To date, 4 studies have tested for an association between the genetic component of gene expression and depression (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ,
      • Gaspar H.A.
      • Gerring Z.
      • Hübel C.
      • Middeldorp C.M.
      • Derks E.M.
      • Breen G.
      Using genetic drug-target networks to develop new drug hypotheses for major depressive disorder.
      ,
      • Gerring Z.F.
      • Gamazon E.R.
      • Derks E.M.
      Major Depressive Disorder Working Group of the Psychiatric Genomics Consortium
      A gene co-expression network-based analysis of multiple brain tissues reveals novel genes and molecular pathways underlying major depression.
      ,
      • Gamazon E.R.
      • Zwinderman A.H.
      • Cox N.J.
      • Denys D.
      • Derks E.M.
      Multi-tissue transcriptome analyses identify genetic mechanisms underlying neuropsychiatric traits.
      ). Wray et al. (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ) showed transcriptomic differences in MD at 17 genes in the dorsolateral prefrontal cortex (DLPFC). Other studies extended these findings by identifying a greater number of associations, for MD (
      • Gaspar H.A.
      • Gerring Z.
      • Hübel C.
      • Middeldorp C.M.
      • Derks E.M.
      • Breen G.
      Using genetic drug-target networks to develop new drug hypotheses for major depressive disorder.
      ,
      • Gerring Z.F.
      • Gamazon E.R.
      • Derks E.M.
      Major Depressive Disorder Working Group of the Psychiatric Genomics Consortium
      A gene co-expression network-based analysis of multiple brain tissues reveals novel genes and molecular pathways underlying major depression.
      ) as well as broad depression (
      • Gamazon E.R.
      • Zwinderman A.H.
      • Cox N.J.
      • Denys D.
      • Derks E.M.
      Multi-tissue transcriptome analyses identify genetic mechanisms underlying neuropsychiatric traits.
      ), across multiple tissues (i.e., distinct brain tissues and whole blood).
      Yet, no TWAS to date has extensively examined these associations. This impeded the identification of which transcriptomic changes are most relevant to depression. As genes correlate in their expression and LD can determine the detection of noncausal genes in TWASs (
      • Mancuso N.
      • Freund M.K.
      • Johnson R.
      • Shi H.
      • Kichaev G.
      • Gusev A.
      • Pasaniuc B.
      Probabilistic fine-mapping of transcriptome-wide association studies.
      ), it is important to better delineate the genetic component of gene expression. Moreover, previous studies generally examined brain and blood tissues only, thus overlooking the role of the hypothalamic-pituitary-adrenal (HPA) and hypothalamic-pituitary-thyroid (HPT) axes. These are important in MD owing to their involvement in the hypo- and hyperactivity of stress responses, sleeping difficulties, and weight loss—all physiological dysfunctions shown in patients with depression (
      • Lopez-Duran N.L.
      • Kovacs M.
      • George C.J.
      Hypothalamic-pituitary-adrenal axis dysregulation in depressed children and adolescents: A meta-analysis.
      ,
      • Stetler C.
      • Miller G.E.
      Depression and hypothalamic-pituitary-adrenal activation: A quantitative summary of four decades of research.
      ,
      • Min W.
      • Liu C.
      • Yang Y.
      • Sun X.
      • Zhang B.
      • Xu L.
      • Sun X.
      Alterations in hypothalamic-pituitary-adrenal/thyroid (HPA/HPT) axes correlated with the clinical manifestations of depression.
      ).
      In the present study, we aimed to delineate key transcriptomic changes within brain, blood, and HPA/HPT axes tissues in depression by exploring the origin of identified transcriptomic associations, in terms of which genes are likely causal and which genes determine both transcriptomic and phenotypic changes (i.e., pleiotropy). Moreover, we examine how such genes lead the dysregulation of nearby coexpressed genes and of biological pathways. This study will add to the field by aiding the understanding of the biological mechanisms of depression, incorporating several disorder-relevant tissues, and highlighting key genes for the disorder that might be useful candidates for future research on its etiology and treatment.

      Methods and Materials

      Datasets

      The analyses used 1) genome-wide summary statistics from the GWAS of MD by Wray et al. (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ), 2) 21 SNP weight sets from 5 separate transcriptomic reference samples, and 3) the 1000 Genomes Project reference for LD estimation.
      First, we leveraged summary statistics from the Psychiatric Genomics Consortium (PGC) MD GWAS (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ), which includes information on the genetic susceptibility to MD for 1,185,038 HapMap3 SNPs from 7 samples: the PGC studies, deCODE, Generation Scotland, GERA (Genetic Epidemiology Research on Adult Health and Aging), iPSYCH, UK Biobank, and 23andMe (not publicly accessible). Participants (n = 135,458 cases, n = 344,901 controls) were of European genetic ancestry.
      Second, SNP weights from distinct tissues and samples (of European genetic ancestry) were used. SNP weights represent the correlations of SNPs with the expression of their annotated gene (
      • Gusev A.
      • Ko A.
      • Shi H.
      • Bhatia G.
      • Chung W.
      • Penninx B.W.J.H.
      • et al.
      Integrative approaches for large-scale transcriptome-wide association studies.
      ). SNP weights from postmortem brain tissue; whole blood; peripheral blood; and adrenal, pituitary, and thyroid glands were downloaded via the TWAS FUSION website (http://gusevlab.org/projects/fusion/#reference-functional-data) after selection based on previous literature (Supplement 1). The weights pertained to the following 5 RNA reference samples: NTR (Netherlands Twin Register) and YFS (Young Finns Study), both of which provide information on blood tissue gene expression; CMC (CommonMind Consortium) and PsychENCODE Consortium, both of which assessed the DLPFC; and the GTEx Consortium, a state-of-the-art study in which expression in multiple brain and peripheral tissues was measured (
      • Gusev A.
      • Ko A.
      • Shi H.
      • Bhatia G.
      • Chung W.
      • Penninx B.W.J.H.
      • et al.
      Integrative approaches for large-scale transcriptome-wide association studies.
      ,
      • Gusev A.
      • Mancuso N.
      • Won H.
      • Kousi M.
      • Finucane H.K.
      • Reshef Y.
      • et al.
      Transcriptome-wide association study of schizophrenia and chromatin activity yields mechanistic disease insights.
      ,
      GTEx Consortium
      Human genomics. The Genotype-Tissue Expression (GTEx) pilot analysis: Multitissue gene regulation in humans.
      ,
      • 1000 Genomes Project Consortium
      • Auton A.
      • Brooks L.D.
      • Durbin R.M.
      • Garrison E.P.
      • Kang 3 H.M.
      • et al.
      A global reference for human genetic variation.
      ), although in a limited number of individuals (determining fewer heritable genes). The SNP weights obtained from a given sample and tissue (e.g., GTEx thyroid, NTR peripheral blood) are called SNP weight sets. Each gene within a given SNP weight set is a feature, i.e., a gene that was examined within a given tissue and sample (e.g., NEGR1 GTEx thyroid). SNP weight set characteristics are presented in Table S1 in Supplement 2.
      Third, the 1000 Genomes Phase 3 European LD reference (N = 489) (
      • Pain O.
      • Pocklington A.J.
      • Holmans P.A.
      • Bray N.J.
      • O’Brien H.E.
      • Hall L.S.
      • et al.
      Novel insight into the etiology of autism spectrum disorder gained by integrating expression data with genome-wide association statistics.
      ) was downloaded from the FUSION website (http://gusevlab.org/projects/fusion/).

      Statistical Analyses

      All statistical analyses were performed in Bash (GNU Project Bourne Again SHell) or R, version 3.5.0 (The R Project for Statistical Computing, Vienna, Austria). Codes and outputs are publicly available at https://opain.github.io/MDD-TWAS.

      Transcriptome-wide Significance Threshold

      To compute the transcriptome-wide significance threshold for this study, we leveraged a permutation procedure used in a previous TWAS (Supplement 1) (
      • Pain O.
      • Pocklington A.J.
      • Holmans P.A.
      • Bray N.J.
      • O’Brien H.E.
      • Hall L.S.
      • et al.
      Novel insight into the etiology of autism spectrum disorder gained by integrating expression data with genome-wide association statistics.
      ). This approach estimates a significance threshold adjusted for the number of tested features, accounting for the correlation between features within and across SNP weight sets. The threshold for transcriptome-wide significance was p = 1.37 × 10−6 (for a false-positive rate of α = .05), with a more stringent significance threshold of p = 3.69 × 10−8, for α = .001.

      TWAS FUSION and Colocalization

      A TWAS FUSION analysis was run on autosomal chromosomes, following the TWAS FUSION protocol with default settings (Supplement 1) (http://gusevlab.org/projects/fusion/). Colocalization was assessed using the coloc R package for all genes meeting transcriptome-wide significance and within a 1.5-Mb window. This Bayesian approach estimates the posterior probability (PP) that associations within a locus for two outcomes are driven by a shared causal variant. It thus enables the distinction between associations driven by horizontal pleiotropy (1 causal SNP affecting both transcription and MD; posterior probability PP4) and linkage (2 causal SNPs in LD affecting transcription and MD separately; posterior probability PP3). More details are available in Supplement 1.

      Conditional Analysis

      A conditional analysis was used to determine whether multiple significant features within a given locus represent independent associations or a single association owing to correlated predicted expression between features. This was performed using the FUSION software, which jointly estimates the effect of all significant features within each locus by using residual SNP associations with depression after accounting for the predicted expression of other features. This process identifies which features represent independent associations (termed jointly significant) and which features are not significant when accounting for the predicted expression of other features in the region (termed marginally significant) (Supplement 1) (
      • Gusev A.
      • Ko A.
      • Shi H.
      • Bhatia G.
      • Chung W.
      • Penninx B.W.J.H.
      • et al.
      Integrative approaches for large-scale transcriptome-wide association studies.
      ). We additionally calculated to what extent the GWAS associations within each locus could be explained by functional associations detected in this TWAS (i.e., the variance explained) (Supplement 1).

      TWAS Fine Mapping

      FOCUS, a TWAS association fine mapping method, was used to identify which features are likely to be causal within regions of association (
      • Mancuso N.
      • Freund M.K.
      • Johnson R.
      • Shi H.
      • Kichaev G.
      • Gusev A.
      • Pasaniuc B.
      Probabilistic fine-mapping of transcriptome-wide association studies.
      ). Analogous to statistical fine mapping of GWAS results, FOCUS estimates the posterior inclusion probability (PIP) of each feature being causal within a region of association, using the sum of PPs to define the default 90% credible set, a set of features likely to contain the causal feature. The method includes a null model where the causal feature is not present. The PIP of individual features is also of interest, with values > .5 indicating that a feature is more likely to be causal than any other feature in the region. The FOCUS fine mapping function was applied across all SNP weight panels simultaneously without the tissue prioritization option.

      TWAS-Based Gene Set Enrichment Analysis

      Finally, we used a previously applied approach for TWAS-based gene set enrichment analysis (TWAS-GSEA) (
      • Pain O.
      • Pocklington A.J.
      • Holmans P.A.
      • Bray N.J.
      • O’Brien H.E.
      • Hall L.S.
      • et al.
      Novel insight into the etiology of autism spectrum disorder gained by integrating expression data with genome-wide association statistics.
      ) to identify functionally informed dysregulated pathways characterizing MD (Supplement 1). A linear mixed model was used to test for an association between z scores indicating nonzero association for each feature and gene set membership, while adjusting for gene length and numbers of SNPs within the gene region and accounting for correlation between features. Linear mixed models were fitted using the R package lme4qtl (
      • Ziyatdinov A.
      • Vázquez-Santiago M.
      • Brunel H.
      • Martinez-Perez A.
      • Aschard H.
      • Soria J.M.
      lme4qtl: Linear mixed models with flexible covariance structure for genetic studies of related individuals.
      ). For the TWAS-GSEA, we used 7246 hypothesis-free gene sets from MSigDB v6.1 (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp) and 76 candidate gene sets from the GWAS by Wray et al. (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ). A minimum of 5 genes within the gene set was required to perform the analysis. TWAS-GSEA was run using different sets of TWAS results to identify gene sets enriched across all tissues, within tissue group (brain, blood, HPA/HPT axes), and within each SNP weight set. If multiple features for a single gene were present, the feature explaining the largest amount of variance in the expression of the gene was retained. Multiple testing correction was applied (false discovery rate q = .05). Moreover, by using BRAINSPAN data, we investigated the preferential expression of genes at multiple developmental stages (19 stages) (
      • Kang H.J.
      • Kawasawa Y.I.
      • Cheng F.
      • Zhu Y.
      • Xu X.
      • Li M.
      • et al.
      Spatio-temporal transcriptome of the human brain.
      ). This was also done within the mixed model method for TWAS-GSEA, in line with previous literature (
      • Pain O.
      • Pocklington A.J.
      • Holmans P.A.
      • Bray N.J.
      • O’Brien H.E.
      • Hall L.S.
      • et al.
      Novel insight into the etiology of autism spectrum disorder gained by integrating expression data with genome-wide association statistics.
      ). A Bonferroni significance threshold was used (p < .002).

      Comparisons With Previous Literature

      Lastly, we compared our findings with previous literature of observed and predicted gene expression. The comparison was performed to evaluate how the findings from the largest study of observed expression in depression to date (N = 1848) (
      • Jansen R.
      • Penninx B.W.J.H.
      • Madar V.
      • Xia K.
      • Milaneschi Y.
      • Hottenga J.J.
      • et al.
      Gene expression in major depressive disorder.
      ) compared with ours. This assessed whole-blood gene expression using microarray technology in the Netherlands Study of Depression and Anxiety. Genes were considered as confirmed if their p value surpassed a Bonferroni-corrected significance threshold accounting for the number of genes compared across studies (p < .05/number of unique genes of nominal significance) in our study. A consistent direction of effect was not a requirement for confirmation of findings. We additionally contrasted our results to 3 previous studies (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ,
      • Gaspar H.A.
      • Gerring Z.
      • Hübel C.
      • Middeldorp C.M.
      • Derks E.M.
      • Breen G.
      Using genetic drug-target networks to develop new drug hypotheses for major depressive disorder.
      ,
      • Gerring Z.F.
      • Gamazon E.R.
      • Derks E.M.
      Major Depressive Disorder Working Group of the Psychiatric Genomics Consortium
      A gene co-expression network-based analysis of multiple brain tissues reveals novel genes and molecular pathways underlying major depression.
      ) of predicted gene expression in MD to evaluate the novel contribution of this study. Overall, we considered 136 unique genes, significant in any SNP weight set, in either of the TWASs.

      Results

      TWAS Analysis

      We identified 176 significant features, from 94 unique genes, which were differentially expressed (p < 1.37 × 10−6) across multiple SNP weight sets (i.e., tissues measured within a sample) in MD (Figure 1; Figure S1A, B in Supplement 1; Table S2 in Supplement 2). Of the 176 significant features, 94 were upregulated, while 82 were downregulated. The most significant feature was NEGR1 (GTEx whole blood) (z = 8.76, p = 1.94 × 10−18). Compared with the GWAS by Wray et al. (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ), 48 unique genes were novel (based on >500 kb distance and R2 < .1) (Table S2 in Supplement 2).
      Figure thumbnail gr1
      Figure 1The relationship between gene expression and major depression. Manhattan-style plot of z scores for each of the tested genes, across all autosomes and tested single nucleotide polymorphism weight sets. Blue lines indicate the transcriptome-wide significance threshold. The names of statistically significant genes are shown.
      The largest number of associations were from the PsychENCODE DLPFC set (22 associated features), but inferences on tissue enrichment are difficult, as SNP weight sets differ in their characteristics (e.g., sample size). Associations are nevertheless shown by tissue group (brain, blood, HPA/HPT axes) in Figure S2A–D in Supplement 1. Of note, 55 associations (from 26 unique genes) were within the extended major histocompatibility complex region (chromosome 6: 26–34 Mb). This region is gene rich and characterized by extensive LD, so these associations should be interpreted with caution. Information on the results once the major histocompatibility complex region is excluded is presented in Supplement 1.

      Colocalization

      We evaluated the colocalization status of a gene by calculating the PP that the genetic and functional associations derived from distinct causal SNPs (PP3) and a shared causal SNP (PP4) (
      • Gusev A.
      • Lawrenson K.
      • Lin X.
      • Lyra P.C.
      • Kar S.
      • Vavra K.C.
      • et al.
      A transcriptome-wide association study of high-grade serous epithelial ovarian cancer identifies new susceptibility genes and splice variants.
      ). Of the 176 significant features, 97 (53 unique genes) were considered as colocalized based on their high PP4 (> .8), in line with previous literature (
      • Taylor K.
      • Davey Smith G.
      • Relton C.L.
      • Gaunt T.R.
      • Richardson T.G.
      Prioritizing putative influential genes in cardiovascular disease susceptibility by applying tissue-specific Mendelian randomization.
      ,
      • Giambartolomei C.
      • Zhenli Liu J.
      • Zhang W.
      • Hauberg M.
      • Shi H.
      • Boocock J.
      • et al.
      A Bayesian framework for multiple trait colocalization from summary association statistics.
      ) (Table S3 in Supplement 2). This means that the same genetic variants were driving associations with depression and with these 97 features, potentially suggesting transcriptomic mediation of genetic risk for depression.

      Conditional Analysis

      We observed that multiple significant features resided within the same locus (defined as a 1.5 Mb ± 0.5 window), for a total of 36 genomic regions (Table 1; Table S4 in Supplement 2). Conditional analysis of the 176 significant features identified 50 jointly (49 unique genes) and 126 marginally (45 unique genes) significant features. This indicated that most of the identified features were associated with depression owing to their coexpression with the 50 independent features. Importantly, anomalous results were shown at 3 loci (including major histocompatibility complex loci), likely reflecting technical problems owing to the mismatch between the LD references used in the different samples (GWAS and functional data) (Table S4 in Supplement 2).
      Table 1Conditional Analysis Findings: Jointly and Marginally Significant Features per Locus
      LocationJointly Significant Features (SNP Weight Set)No. of Marginally Significant FeaturesTop TWAS pTop GWAS pVariance Explained, %
      Variance in the top GWAS SNP explained by the top TWAS feature within the locus.
      chr1:7413452-9875347RERE (YFS blood)31.10 × 10−73.18 × 10−8100.00
      chr1:35885799-37876701SNORA63 (GTEx nucleus accumbens)01.24 × 10−66.27 × 10−869.30
      chr1:71753372-73766162NEGR1 (GTEx whole blood)51.94 × 10−184.54 × 10−1597.80
      chr1:174891875-177103690RFWD2 (CMC DLPFC splicing)34.66 × 10−72.30 × 10−790.00
      chr1:180725304-182724241CACNA1E (CMC DLPFC splicing)06.06 × 10−71.08 × 10−764.00
      chr1:196478918-198741422DENND1B (CMC DLPFC splicing)25.90 × 10−83.11 × 10−892.60
      chr2:57388379-59467945FANCL (CMC DLPFC splicing, CMC DLPFC)02.18 × 10−74.68 × 10−985.00
      chr2:196832647-199295649ANKRD44 (YFS blood)11.27 × 10−83.52 × 10−782.80
      chr3:43487406-45561063ZNF445 (CMC DLPFC)03.34 × 10−76.34 × 10−874.70
      chr4:40937584-43090938SLC30A9 (GTEx cortex), TMEM33 (PsychENCODE DLPFC)87.72 × 10−93.59 × 10−992.30
      chr5:139030460-141219083PCDHA5 (GTEx thyroid)26.55 × 10−81.37 × 10−687.70
      chr6:98832858-100829135COQ3 (CMC DLPFC splicing)02.65 × 10−89.09 × 10−835.10
      chr6:104405706-106583999BVES-AS1 (GTEx amygdala)22.43 × 10−89.50 × 10−892.90
      chr7:11252396-13282905TMEM106B (PsychENCODE DLPFC)37.01 × 10−92.55 × 10−8100.00
      chr7:24021857-26019767OSBPL3 (GTEx pituitary)01.88 × 10−86.49 × 10−777.70
      chr8:51238261-53720740PXDNL (CMC DLPFC)03.92 × 10−91.34 × 10−783.80
      chr8:60435234-62428932RP11-163N6.2 (GTEx thyroid)09.47 × 10−85.25 × 10−789.80
      chr9:125606617-127604411PIGFP2 (PsychENCODE DLPFC)01.12 × 10−72.73 × 10−863.80
      chr11:56092913-58422547TNKS1BP1 (GTEx adrenal gland), CLP1 (GTEx whole blood)12.04 × 10−71.47 × 10−795.20
      chr11:60540194-62557903TMEM258 (PsychENCODE DLPFC)05.12 × 10−74.26 × 10−783.90
      chr11:112346414-114345882DRD2 (GTEx frontal cortex)03.90 × 10−74.90 × 10−70.41
      chr13:52652520-54625616OLFM4 (CMC DLPFC)03.56 × 10−76.06 × 10−1929.90
      chr14:41077086-43073683CTD-2298J14.2 (GTEx thyroid)21.36 × 10−82.57 × 10−988.10
      chr14:58952573-61334943CCDC175 (GTEx thyroid)34.28 × 10−82.18 × 10−782.30
      chr14:63322572-65770213ESR2 (GTEx pituitary)22.20 × 10−97.60 × 10−1080.00
      chr14:74120633-76388050PROX2 (GTEx thyroid)58.51 × 10−96.71 × 10−993.50
      chr14:102878783-105180229RP11-894P9.2 (GTEx thyroid)114.69 × 10−83.05 × 10−984.60
      chr16:71147494-73210261PMFBP1 (PsychENCODE DLPFC)02.46 × 10−73.35 × 10−876.30
      chr17:26406423-28478661TIAF1 (GTEx adrenal gland)18.27 × 10−88.51 × 10−958.50
      chr17:64524284-66521332CTD-2653B5.1 (PsychENCODE DLPFC)03.30 × 10−75.39 × 10−625.80
      chr18:51385406-53561919RAB27B (PsychENCODE DLPFC)15.36 × 10−73.62 × 10−1114.60
      chr20:46838019-48853908DDX27 (CMC DLPFC)01.32 × 10−63.54 × 10−691.00
      chr22:40218102-42697216ZC3H7B (GTEx cerebellum)81.01 × 10−87.56 × 10−995.50
      Chr, chromosome; CMC, CommonMind Consortium; DLPFC, dorsolateral prefrontal cortex; GTEx, Genotype-Tissue Expression; GWAS, genome-wide association study; NTR, Netherlands Twin Register; SNP, single nucleotide polymorphism; TWAS, transcriptome-wide association study; YFS, Young Finns Study.
      a Variance in the top GWAS SNP explained by the top TWAS feature within the locus.
      We additionally investigated the effect of adjusting for features’ gene expression on associations between SNPs and the trait (i.e., GWAS findings). The variance in the GWAS associations accounted for by gene expression associations at a given locus ranged from 0.41% to 100%, with a median variance explained of 84%. Genome-wide associations at 3 loci, including TMEM106B (PsychENCODE DLPFC), were fully explained by our TWAS results (Table 1). This might mean that at these loci, genetic risk for depression is mediated by functional changes.

      Statistical Fine Mapping

      FOCUS was performed to calculate the probability estimates of causality (PIP) for each feature. We found 23 features (23 unique genes) with PIP > .5, indicating these are likely causal in their associations with depression (Table 2; Table S5 in Supplement 2). Of these, 11 were supported by the colocalization analysis (Table 2). The highest probability of causality was within NEGR1 GTEx whole blood (PIP = 1.00).
      Table 2Statistical Fine Mapping Findings: Potentially Causal Features
      LocationGeneSNP Weight SetFOCUS PIP
      PIP estimate of causality.
      Colocalization
      chr1:8412457-8877702REREYFS blood.50True
      chr1:36884051-36884179SNORA63GTEx nucleus accumbens.70False
      chr1:71861623-72748417NEGR1GTEx whole blood1.00True
      chr1:181452685-181775921CACNA1ECMC DLPFC splicing.52False
      chr1:197473878-197744623DENND1BCMC DLPFC.73True
      chr2:58386377-58468515FANCLCMC DLPFC splicing.80True
      chr4:41992489-42092474SLC30A9GTEx cortex.53False
      chr5:87988462-87989789CTC-467M3.3GTEx frontal cortex.84True
      chr6:99817347-99842082COQ3CMC DLPFC splicing.97False
      chr6:105584224-105617820BVES-AS1GTEx amygdala.74False
      chr7:12250867-12282993TMEM106BGTEx whole blood.61True
      chr7:24836158-25021253OSBPL3GTEx pituitary.98False
      chr8:52232136-52722005PXDNLCMC DLPFC1.00False
      chr8:61297147-61429354RP11-163N6.2GTEx thyroid.92False
      chr9:126605315-126605965PIGFP2PsychENCODE DLPFC.94False
      chr11:113280318-113346111DRD2GTEx frontal cortex.97False
      chr13:53602875-53626196OLFM4CMC DLPFC.99False
      chr14:42076773-42373752LRFN5GTEx cerebellum.50True
      chr14:59971779-60043549CCDC175GTEx thyroid.61True
      chr14:64550950-64770377ESR2GTEx pituitary.59True
      chr14:75319736-75330537PROX2GTEx thyroid.72True
      chr16:72146056-72210777PMFBP1PsychENCODE DLPFC.96False
      chr17:27401933-27405875TIAF1GTEx adrenal gland.75True
      Chr, chromosome; CMC, CommonMind Consortium; DLPFC, dorsolateral prefrontal cortex; GTEx, Genotype-Tissue Expression; PIP, posterior inclusion probability; SNP, single nucleotide polymorphism; YFS, Young Finns Study.
      a PIP estimate of causality.

      High-Confidence Associations

      In a data-driven approach, we highlighted the genes that are most relevant to MD (high-confidence associations). First, in the TWAS results, we applied a more stringent significance threshold of p < 3.69 × 10−8 (α = .001) to minimize the chance of false-positive results (
      • Colquhoun D.
      An investigation of the false discovery rate and the misinterpretation of p-values.
      ). Second, from these highly significant genes, we identified the genes that were colocalized (PP4 > .8) and likely to be causal (PIP > .5). These features were therefore strongly associated with MD, possibly at a causal level, and their transcription changes resulted from the corresponding genetic risk for depression. This strategy showed 6 high-confidence associations (Table 3), from brain tissues (frontal cortex, cerebellum) and peripheral tissues (whole blood, thyroid, pituitary). The NEGR1 gene (GTEx whole blood) constituted the most significant hit (p = 1.94 × 10−18) with the highest probability of causality (PIP = 1.00).
      Table 3Characteristics of High-Confidence Associations: Highly Significant, Colocalized, and Potentially Causal Features
      LocationGeneSNP Weight SetTWAS zTWAS pNovel
      Genes are novel compared with GWAS if >500 kb away from a lead GWAS variant and if their predicted expression is not correlated with GWAS variant (R2 < .1). Genes are novel compared with previous TWASs of major depression if these did not identify them as statistically significant.
      Colocalization PP4
      PP of colocalization.
      FOCUS PIP
      PIP estimate of causality.
      chr1:71861623-72748417NEGR1GTEx whole blood8.761.94 × 10−18No.991.00
      chr5:87988462-87989789CTC-467M3.3GTEx frontal cortex−7.091.33 × 10−12No.85.84
      chr7:12250867-12282993TMEM106BGTEx whole blood5.533.18 × 10−8Yes (GWAS).99.61
      chr14:42076773-42373752LRFN5GTEx cerebellum5.602.17 × 10−8No.96.50
      chr14:64550950-64770377ESR2GTEx pituitary−5.982.20 × 10−9No.86.59
      chr14:75319736-75330537PROX2GTEx thyroid−5.768.51 × 10−9Yes (TWAS).96.72
      Chr, chromosome; GTEx, Genotype-tissue expression; GWAS, genome-wide association study; PIP, posterior inclusion probability; PP, posterior probability; SNP, single nucleotide polymorphism; TWAS, transcriptome-wide association study.
      a Genes are novel compared with GWAS if >500 kb away from a lead GWAS variant and if their predicted expression is not correlated with GWAS variant (R2 < .1). Genes are novel compared with previous TWASs of major depression if these did not identify them as statistically significant.
      b PP of colocalization.
      c PIP estimate of causality.

      TWAS-Based Gene Set Enrichment Analysis

      Candidate gene analysis highlighted enrichment for 15 gene sets, including neuronal and synaptic processes, schizophrenia genetic risk, and RBFOX2, a key regulator of splicing within the nervous system and of estrogen receptor transcriptional activity (Table 4) (
      • Geer L.Y.
      • Marchler-Bauer A.
      • Geer R.C.
      • Han L.
      • He J.
      • He S.
      • et al.
      The NCBI BioSystems database.
      ). The hypothesis-free analysis revealed enrichment for 7 gene sets, including macromolecular and protein complex binding (Table 5). Enrichment or depletion of our differentially expressed genes was not found at any developmental stage (Figure S3A–C in Supplement 1).
      Table 4Candidate Gene Set Enrichment Results Based on TWAS Results From All Tissues, Tissue Sets, and Each Panel
      Gene SetPubMed IDEstimateSEtpFDR p
      Multiple testing–corrected p value using FDR method.
      Tissue
      All Tissue Results
      RBFOX2246133500.080.023.265.57 × 10−43.95 × 10−2All
      Tissue-Set Results
      RBFOX2246133500.110.034.131.78 × 10−51.26 × 10−3Brain
      SCZ.COMPOSITE244635080.110.033.641.39 × 10−44.92 × 10−3Brain
      RBFOX1.RBFOX3246133500.080.023.452.80 × 10−46.62 × 10−3Brain
      FMRP217842460.120.043.353.98 × 10−47.06 × 10−3Brain
      POTENTIALLY.SYNAPTIC.ALL276949940.060.023.031.21 × 10−31.72 × 10−2Brain
      PGC.BP.P10.4219269720.180.062.852.20 × 10−32.60 × 10−2Brain
      NEURONAL.PSD230716130.090.032.673.83 × 10−33.89 × 10−2Brain
      Tissue-Specific Results
      MIR.137244635080.360.103.581.70 × 10−41.02 × 10−2CMC DLPFC RNA-seq
      SCZ.DENOVO.NONSYN244635080.400.123.512.25 × 10−41.17 × 10−2GTEx pituitary
      SCZ.COMPOSITE244635080.230.073.285.21 × 10−41.36 × 10−2GTEx pituitary
      SCZ.COMPOSITE244635080.270.083.197.03 × 10−43.24 × 10−2GTEx basal ganglia
      CONSTRAINED250866660.270.092.971.50 × 10−33.44 × 10−2CMC DLPFC RNA-seq
      RBFOX1.RBFOX3246133500.120.042.931.72 × 10−33.44 × 10−2CMC DLPFC RNA-seq
      PGC.SCZ.P10.4244635080.270.102.733.20 × 10−34.79 × 10−2CMC DLPFC RNA-seq
      CMC, CommonMind Consortium; DLPFC, dorsolateral prefrontal cortex; FDR, false discovery rate; GTEx, Genotype Tissue Expression; RNA-seq, RNA sequencing; SE, standard error; TWAS, transcriptome-wide association study.
      a Multiple testing–corrected p value using FDR method.
      Table 5Hypothesis-free Gene Set Enrichment Results Based on TWAS Results From Each Panel
      Gene SetEstimateSEtpFDR p
      Multiple testing–corrected p value using FDR method.
      Tissue
      GO.MACROMOLECULAR.COMPLEX.BINDING0.460.095.071.98 × 10−75.99 × 10−4GTEx basal ganglia
      GO.MICROTUBULE.BINDING1.040.224.681.46 × 10−62.20 × 10−3GTEx basal ganglia
      GO.ALCOHOL.BINDING1.790.384.661.56 × 10−65.16 × 10−3GTEx pituitary
      GO.CHROMATIN.BINDING0.790.194.191.42 × 10−51.41 × 10−2GTEx basal ganglia
      GO.PROTEIN.COMPLEX.BINDING0.440.114.121.87 × 10−51.41 × 10−2GTEx basal ganglia
      GO.LIGAND.DEPENDENT.NUCLEAR.RECEPTOR.BINDING1.980.503.983.47 × 10−52.09 × 10−2GTEx basal ganglia
      GO.REGULATION.OF.INTRINSIC.APOPTOTIC.SIGNALING.PATHWAY1.790.434.141.73 × 10−53.30 × 10−2GTEx amygdala
      FDR, false discovery rate; GTEx, Genotype Tissue Expression; SE, standard error; TWAS, transcriptome-wide association study.
      a Multiple testing–corrected p value using FDR method.

      Comparison With Previous Literature

      First, our findings were compared with the largest study of observed gene expression to date for MD (
      • Jansen R.
      • Penninx B.W.J.H.
      • Madar V.
      • Xia K.
      • Milaneschi Y.
      • Hottenga J.J.
      • et al.
      Gene expression in major depressive disorder.
      ). Of the significant genes (false discovery rate < 0.1) of such study, we also identified 6 unique genes (PAPPA2, MBNL1, TMEM64, GNPTAB, KTN1, TMED10). Four of these (PAPPA2, GNPTAB, KTN1, TMED10) were consistently upregulated or downregulated (Table S6 in Supplement 1; Supplement 1).
      To evaluate the novelty of our findings, we compared our results with previous studies of predicted gene expression changes in MD performed by Wray et al. (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ) (TWAS FUSION method), Gaspar et al. (
      • Gaspar H.A.
      • Gerring Z.
      • Hübel C.
      • Middeldorp C.M.
      • Derks E.M.
      • Breen G.
      Using genetic drug-target networks to develop new drug hypotheses for major depressive disorder.
      ) (S-PrediXcan), and Gerring et al. (
      • Gerring Z.F.
      • Gamazon E.R.
      • Derks E.M.
      Major Depressive Disorder Working Group of the Psychiatric Genomics Consortium
      A gene co-expression network-based analysis of multiple brain tissues reveals novel genes and molecular pathways underlying major depression.
      ) (S-PrediXcan). Our results were overlapping with previous TWASs at 83 genes (Table S7 in Supplement 2). Moreover, across all TWASs performed to date, we identified 53 novel associations (unique genes). This might be due to differences in the type and number of SNP weight sets used, methods, and statistical thresholds (e.g., permutation-based vs. Bonferroni).

      Discussion

      This is the first study to uncover the genetic component of gene expression in MD through a comprehensive investigation in multiple tissues and an in-depth characterization of the identified associations delineating key transcriptomic changes. Here, we highlight a few key findings. First, we detected 94 unique genes associated with depression, approximately half of which were novel. Second, with a series of follow-up analyses, we found 6 associations that were of high confidence, based on their strong TWAS associations, colocalization, and probability of causality. These findings highlight the transcriptomic changes in MD that likely play a role in etiology of the disorder, resulting from the same genetic variations as the disorder and determining, in turn, transcriptomic changes at nearby genes and pathways involved in key biological processes. Such findings on key genes for depression can guide future research on drug targets as well as candidate gene investigations in animal studies, where consequences of molecular alterations can be more readily observed by inducing gene knockdown or upregulation.

      Key Findings

      When testing for an association between gene expression and MD, we detected 94 transcriptome-wide significant genes, differentially expressed across multiple tissues, thus demonstrating the presence of widespread transcriptomic changes in depression. Comparison with previous literature highlighted the novelty of this study, which enabled the identification of 48 (compared with GWAS findings) and 53 (compared with previous TWASs) novel genes.
      Further investigation of significant associations through a conditional analysis determined whether gene associations within the same genomic region represent independent associations or whether multiple genes are associated owing to correlated predicted expression. The 94 significant genes represented 49 independent associations with depression, suggesting that approximately half of the identified signal depends on LD and correlated predicted expression of nearby genes. The strength of association for each feature can be affected by different characteristics of the SNP weight sets (e.g., sample size), thus warranting cautionary interpretations of which features are driving the associated loci. Comparison of the GWAS summary statistics before and after conditioning on significant TWAS associations additionally revealed that GWAS associations were explained to a major extent by TWAS associations, further suggesting the possibility of transcriptomic mediation of genetic risk for depression.
      When exploring whether significant associations were driven by pleiotropy or linkage using a colocalization analysis, we observed that 53 transcription-MD relationships derived from the same causal polymorphisms underlying SNP-MD associations. This indicated that most of the detected genes constituted pleiotropic effects as opposed to linkage. While these findings suggest that transcription mediates the relationship between genetic susceptibility and depression, colocalization does not test for such relationships, and it cannot identify specific causal variants (
      • Gusev A.
      • Ko A.
      • Shi H.
      • Bhatia G.
      • Chung W.
      • Penninx B.W.J.H.
      • et al.
      Integrative approaches for large-scale transcriptome-wide association studies.
      ).
      To gain further insight into which genes are likely causal for MD, we used a TWAS fine mapping approach called FOCUS. In some instances, FOCUS clearly highlighted a single feature as the causal association, such as the upregulation of NEGR1 in the blood. Of the 23 genes with a high probability of causality estimated by FOCUS, 11 were identified as colocalized. A key distinction between these methods is that colocalization assumes a single causal variant, whereas FOCUS allows for multiple causal variants.
      Six high-confidence associations were identified. Notably, none of these were found by previous TWASs on other psychiatric phenotypes (e.g., bipolar disorder) (
      • Gamazon E.R.
      • Zwinderman A.H.
      • Cox N.J.
      • Denys D.
      • Derks E.M.
      Multi-tissue transcriptome analyses identify genetic mechanisms underlying neuropsychiatric traits.
      ,
      • Liao C.
      • Laporte A.D.
      • Spiegelman D.
      • Akçimen F.
      • Joober R.
      • Dion P.A.
      • Rouleau G.A.
      Transcriptome-wide association study of attention deficit hyperactivity disorder identifies 6 associated genes and phenotypes.
      ,
      • Huckins L.M.
      • Dobbyn A.
      • Ruderfer D.M.
      • Hoffman G.
      • Wang W.
      • Pardiñas A.F.
      • et al.
      Gene expression imputation across multiple brain regions provides insights into 9 schizophrenia risk.
      ). Of the high-confidence associations, 3 should be highlighted owing to their functional role: NEGR1, ESR2, and TMEM106B. SNPs within these 3 genes have been previously detected in one or more GWASs of depression (
      • Howard D.M.
      • Adams M.J.
      • Clarke T.K.
      • Hafferty J.D.
      • Gibson J.
      • Shirali M.
      • et al.
      Genome-wide meta-analysis of depression identifies 102 independent variants and highlights the importance of the prefrontal brain regions.
      ,
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ,
      • Hyde C.L.
      • Nagle M.W.
      • Tian C.
      • Chen X.
      • Paciga S.A.
      • Wendland J.R.
      • et al.
      Identification of 15 genetic loci associated with risk of major depression in individuals of European descent.
      ). Our study contributed to previous literature by elucidating the functional characteristics of such genes, showing upregulation for NEGR1 and TMEM106B and downregulation for ESR2.
      NEGR1 plays a role in axon extension, synaptic plasticity, and synapse formation, processes key to neuronal functioning (
      • Hashimoto T.
      • Maekawa S.
      • Miyata S.
      IgLON cell adhesion molecules regulate synaptogenesis in hippocampal neurons.
      ,
      • Pischedda F.
      • Szczurkowska J.
      • Cirnaru M.D.
      • Giesert F.
      • Vezzoli E.
      • Ueffing M.
      • et al.
      A cell surface biotinylation assay to reveal membrane-associated neuronal cues: Negr1 regulates dendritic arborization.
      ,
      • Sanz R.
      • Ferraro G.B.
      • Fournier A.E.
      IgLON cell adhesion molecules are shed from the cell surface of cortical neurons to promote neuronal growth.
      ). Moreover, it includes variants found to be related to obesity, a trait repeatedly correlated with MD, at both phenotypic and biological levels (
      • Milaneschi Y.
      • Simmons W.K.
      • Rossum EFC van
      • Penninx B.W.
      Depression and obesity: Evidence of shared biological mechanisms.
      ). ESR2 regulates the activity of estrogen, a sex hormone involved in HPA axis activity and inflammation (
      • Ochedalski T.
      • Subburaju S.
      • Wynn P.C.
      • Aguilera G.
      Interaction between oestrogen and oxytocin on hypothalamic-pituitary-adrenal axis activity.
      ,
      • Puder J.J.
      • Freda P.U.
      • Goland R.S.
      • Wardlaw S.L.
      Estrogen modulates the hypothalamic-pituitary-adrenal and inflammatory cytokine responses to endotoxin in women.
      ). Estrogen fluctuations across a woman’s life span, for example, during premenstrual monthly period and menopause, have been proposed as a risk factor for depression (
      • Soares C.N.
      • Zitek B.
      Reproductive hormone sensitivity and risk for depression across the female life cycle: A continuum of vulnerability?.
      ). However, it remains unclear whether the risk is driven by estrogen or other co-occurring hormonal changes (
      • Soares C.N.
      • Zitek B.
      Reproductive hormone sensitivity and risk for depression across the female life cycle: A continuum of vulnerability?.
      ,
      • Behl C.
      Oestrogen as a neuroprotective hormone.
      ). TMEM106B was previously implicated as a susceptibility gene for neurodegenerative disorders (
      • Hernández I.
      • Rosende-Roca M.
      • Alegret M.
      • Mauleón A.
      • Espinosa A.
      • Vargas L.
      • et al.
      Association of TMEM106B rs1990622 marker and frontotemporal dementia: Evidence for a recessive effect and meta-analysis.
      ,
      • Nicholson A.M.
      • Rademakers R.
      What we know about TMEM106B in neurodegeneration.
      ,
      • Rutherford N.J.
      • Carrasquillo M.M.
      • Li M.
      • Bisceglio G.
      • Menke J.
      • Josephs K.A.
      • et al.
      TMEM106B risk variant is implicated in the pathologic presentation of Alzheimer disease.
      ) and TDP-43 abnormalities (
      • Nicholson A.M.
      • Rademakers R.
      What we know about TMEM106B in neurodegeneration.
      ), which are featured in such pathologies (
      • Nicholson A.M.
      • Rademakers R.
      What we know about TMEM106B in neurodegeneration.
      ,
      • Rutherford N.J.
      • Carrasquillo M.M.
      • Li M.
      • Bisceglio G.
      • Menke J.
      • Josephs K.A.
      • et al.
      TMEM106B risk variant is implicated in the pathologic presentation of Alzheimer disease.
      ,
      • Van Deerlin V.M.
      • Sleiman P.M.A.
      • Martinez-Lage M.
      • Chen-Plotkin A.
      • Wang L.-S.
      • Graff-Radford N.R.
      • et al.
      Common variants at 7p21 are associated with frontotemporal lobar degeneration with TDP-43 inclusions.
      ). Depression has been repeatedly suggested as a risk factor for neurodegenerative disorders (
      • Byers A.L.
      • Yaffe K.
      Depression and risk of developing dementia.
      ,
      • Li J.Q.
      • Tan L.
      • Wang H.F.
      • Tan M.S.
      • Tan L.
      • Xu W.
      • et al.
      Risk factors for predicting progression from mild cognitive impairment to Alzheimer’s disease: A systematic review and meta-analysis of cohort studies.
      ,
      • Ownby R.L.
      • Crocco E.
      • Acevedo A.
      • John V.
      • Loewenstein D.
      Depression and risk for Alzheimer disease: Systematic review, meta-analysis, and metaregression analysis.
      ). Moreover, TDP-43 proteinopathy was shown in a small sample of patients with late-life depression (
      • Ichikawa T.
      • Baba H.
      • Maeshima H.
      • Shimano T.
      • Inoue M.
      • Ishiguro M.
      • et al.
      Serum levels of TDP-43 in late-life patients with depressive episode.
      ). Overall, while previous literature points to an important role of these genes in depression and related phenotypes, replication of associations is necessary as well as greater insight into the relationship between estrogen levels and TDP-43 with MD.
      TWAS-GSEA was able to identify several gene sets showing dysregulated expression in individuals with MD. The enriched gene sets are congruent with a previous GWAS-based enrichment analysis (
      • Wray N.R.
      • Ripke S.
      • Mattheisen M.
      • Trzaskowski M.
      • Byrne E.M.
      • Abdellaoui A.
      • et al.
      Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
      ), corroborating the importance, among others, of genes bound by transcription factors (RBFOX1, RBFOX2, RBFOX3, FMRP) and genes encoding synaptic proteins and ion channels. Novel enriched pathways were also found for key biological functions such as protein and macromolecular complex binding. Replication is warranted.

      Limitations of Present Study and Suggestions for Future Research

      While these findings are promising, several limitations merit discussion. First, the small sample sizes of the gene expression reference samples may have impeded the detection of subtle effects of the transcriptome on depression, meaning that larger samples are needed. Second, the wide range of MD definitions within the GWAS samples, ranging from self-reported depression (e.g., 23andMe) to clinical diagnosis (e.g., iPSYCH), might impact results. Future studies could investigate to what extent minimal phenotyping affects findings at the transcriptomic level compared with more robust definitions. Third, we analyzed a wider set of tissues than previous MD TWASs, with 21 distinct SNP weight sets from blood, brain, and HPA/HPT axes. This may uncover more true associations but may also have introduced noise, as using tissues not strictly relevant to depression might capture noncausal genes (
      • Wainberg M.
      • Sinnott-Armstrong N.
      • Mancuso N.
      • Barbeira A.N.
      • Knowles D.A.
      • Golan D.
      • et al.
      Opportunities and challenges for transcriptome-wide association studies.
      ). Nonetheless, tested tissues were selected based on previous literature, meaning that such tissues are supposedly disease relevant. An alternative strategy would be to use data-driven tissue selection, leveraging an LD score regression-based method (
      • Finucane H.K.
      • Reshef Y.A.
      • Anttila V.
      • Slowikowski K.
      • Gusev A.
      • Byrnes A.
      • et al.
      Heritability enrichment of specifically expressed genes identifies disease-relevant tissues and cell types.
      ). Moreover, valuable eQTL data (e.g., eQTLGen Consortium) were not used because of the lack of precomputed SNP weight and access to individual-level data (for weights development). We nevertheless recommend the calculation of such weights when possible. Furthermore, our TWAS approach, by solely assessing the cis-genetic component of gene expression, cannot capture trans-eQTL effects. Future research should channel resources toward building larger gene expression reference panels to enable investigation of trans-eQTL effects. Lastly, while this study provided further insight into the relationship between SNPs, gene expression, and depression and used colocalization and causal fine mapping analyses to test certain criteria of a causal model, it cannot verify causality between associated genes and depression.
      In conclusion, we provide evidence for widespread transcriptomic changes in MD. Our study enables the detection of novel associations and the elucidation of the transcriptomic changes that previously identified risk genes undergo. We underline genes that might be of key relevance to depression, including NEGR1, ESR2, and TMEM106B. These results suggest an important role of the genetic component of gene expression in depression.

      Acknowledgments and Disclosures

      This work was supported by the Medical Research Council (Grant No. MR/N015746/1 [to CML] ) and National Institute for Health Research Biomedical Research Centre at South London and Maudsley National Health Service Foundation Trust and King’s College London. The authors acknowledge the use of the research computing facility at King’s College London, Rosalind (https://rosalind.kcl.ac.uk), which is delivered in partnership with the National Institute for Health Research Maudsley Biomedical Research Centre and partly funded by capital equipment grants from the Maudsley Charity (Grant No. 980 ) and Guy’s and St. Thomas’ Charity (Grant No. TR130505 ). The views expressed are those of the authors and not necessarily those of the National Health Service, National Institute for Health Research, or Department of Health and Social Care.
      We thank Alexander Gusev for assisting us and for creating the TWAS FUSION method used here. We also thank the research participants and employees of 23andMe for making this work possible.
      A previous version of this article was published as a preprint on bioRxiv: https://doi.org/10.1101/2020.03.24.004903.
      CML sits on the Myriad Neuroscience Scientific Advisory Board. The other authors report no biomedical financial interests or potential conflicts of interest.

      Supplementary Material

      References

        • World Health Organization
        Depression and Other Common Mental Disorders: Global Health Estimates..
        (Available at:) (Accessed November 9, 2020)
        • American Psychiatric Association
        Diagnostic and Statistical Manual of Mental Disorders.
        5th ed. American Psychiatric Association, Arlington, VA2013
        • Sullivan P.F.
        • Neale M.C.
        • Kendler K.S.
        Genetic epidemiology of major depression: Review and meta-analysis.
        Am J Psychiatry. 2000; 157: 1552-1562
        • Howard D.M.
        • Adams M.J.
        • Clarke T.K.
        • Hafferty J.D.
        • Gibson J.
        • Shirali M.
        • et al.
        Genome-wide meta-analysis of depression identifies 102 independent variants and highlights the importance of the prefrontal brain regions.
        Nat Neurosci. 2019; 22: 343
        • Wray N.R.
        • Ripke S.
        • Mattheisen M.
        • Trzaskowski M.
        • Byrne E.M.
        • Abdellaoui A.
        • et al.
        Genome-wide association analyses identify 44 risk variants and refine the genetic architecture of major depression.
        Nat Genet. 2018; 50: 668
        • GTEx Consortium
        Genetic effects on gene expression across human tissues.
        Nature. 2017; 550: 204-213
        • Barbeira A.N.
        • Dickinson S.P.
        • Bonazzola R.
        • Zheng J.
        • Wheeler H.E.
        • Torres J.M.
        • et al.
        Exploring the phenotypic consequences of tissue specific gene expression variation inferred from GWAS summary statistics.
        Nat Commun. 2018; 9: 1825
        • Gusev A.
        • Ko A.
        • Shi H.
        • Bhatia G.
        • Chung W.
        • Penninx B.W.J.H.
        • et al.
        Integrative approaches for large-scale transcriptome-wide association studies.
        Nat Genet. 2016; 48: 245-252
        • Gaspar H.A.
        • Gerring Z.
        • Hübel C.
        • Middeldorp C.M.
        • Derks E.M.
        • Breen G.
        Using genetic drug-target networks to develop new drug hypotheses for major depressive disorder.
        Transl Psychiatry. 2019; 9: 117
        • Gerring Z.F.
        • Gamazon E.R.
        • Derks E.M.
        • Major Depressive Disorder Working Group of the Psychiatric Genomics Consortium
        A gene co-expression network-based analysis of multiple brain tissues reveals novel genes and molecular pathways underlying major depression.
        PLoS Genet. 2019; 15: e1008245
        • Gamazon E.R.
        • Zwinderman A.H.
        • Cox N.J.
        • Denys D.
        • Derks E.M.
        Multi-tissue transcriptome analyses identify genetic mechanisms underlying neuropsychiatric traits.
        Nat Genet. 2019; 51: 933
        • Mancuso N.
        • Freund M.K.
        • Johnson R.
        • Shi H.
        • Kichaev G.
        • Gusev A.
        • Pasaniuc B.
        Probabilistic fine-mapping of transcriptome-wide association studies.
        Nat Genet. 2019; 51: 675-682
        • Lopez-Duran N.L.
        • Kovacs M.
        • George C.J.
        Hypothalamic-pituitary-adrenal axis dysregulation in depressed children and adolescents: A meta-analysis.
        Psychoneuroendocrinology. 2009; 34: 1272-1283
        • Stetler C.
        • Miller G.E.
        Depression and hypothalamic-pituitary-adrenal activation: A quantitative summary of four decades of research.
        Psychosom Med. 2011; 73: 114
        • Min W.
        • Liu C.
        • Yang Y.
        • Sun X.
        • Zhang B.
        • Xu L.
        • Sun X.
        Alterations in hypothalamic-pituitary-adrenal/thyroid (HPA/HPT) axes correlated with the clinical manifestations of depression.
        Prog Neuropsychopharmacol Biol Psychiatry. 2012; 39: 206-211
        • Gusev A.
        • Mancuso N.
        • Won H.
        • Kousi M.
        • Finucane H.K.
        • Reshef Y.
        • et al.
        Transcriptome-wide association study of schizophrenia and chromatin activity yields mechanistic disease insights.
        Nat Genet. 2018; 50: 538-548
        • GTEx Consortium
        Human genomics. The Genotype-Tissue Expression (GTEx) pilot analysis: Multitissue gene regulation in humans.
        Science. 2015; 348: 648-660
        • 1000 Genomes Project Consortium
        • Auton A.
        • Brooks L.D.
        • Durbin R.M.
        • Garrison E.P.
        • Kang 3 H.M.
        • et al.
        A global reference for human genetic variation.
        Nature. 2015; 52668–74
        • Pain O.
        • Pocklington A.J.
        • Holmans P.A.
        • Bray N.J.
        • O’Brien H.E.
        • Hall L.S.
        • et al.
        Novel insight into the etiology of autism spectrum disorder gained by integrating expression data with genome-wide association statistics.
        Biol Psychiatry. 2019; 86: 265-273
        • Ziyatdinov A.
        • Vázquez-Santiago M.
        • Brunel H.
        • Martinez-Perez A.
        • Aschard H.
        • Soria J.M.
        lme4qtl: Linear mixed models with flexible covariance structure for genetic studies of related individuals.
        BMC Bioinformatics. 2018; 19: 68
        • Kang H.J.
        • Kawasawa Y.I.
        • Cheng F.
        • Zhu Y.
        • Xu X.
        • Li M.
        • et al.
        Spatio-temporal transcriptome of the human brain.
        Nature. 2011; 478: 483-489
        • 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: 339-347
        • Gusev A.
        • Lawrenson K.
        • Lin X.
        • Lyra P.C.
        • Kar S.
        • Vavra K.C.
        • et al.
        A transcriptome-wide association study of high-grade serous epithelial ovarian cancer identifies new susceptibility genes and splice variants.
        Nat Genet. 2019; 51: 815-823
        • Taylor K.
        • Davey Smith G.
        • Relton C.L.
        • Gaunt T.R.
        • Richardson T.G.
        Prioritizing putative influential genes in cardiovascular disease susceptibility by applying tissue-specific Mendelian randomization.
        Genome Med. 2019; 11: 6
        • Giambartolomei C.
        • Zhenli Liu J.
        • Zhang W.
        • Hauberg M.
        • Shi H.
        • Boocock J.
        • et al.
        A Bayesian framework for multiple trait colocalization from summary association statistics.
        Bioinformatics. 2018; 34: 2538-2545
        • Colquhoun D.
        An investigation of the false discovery rate and the misinterpretation of p-values.
        R Soc Open Sci. 2014; 1: 140216
        • Geer L.Y.
        • Marchler-Bauer A.
        • Geer R.C.
        • Han L.
        • He J.
        • He S.
        • et al.
        The NCBI BioSystems database.
        Nucleid Acids Res. 2010; 38: D492-D496
        • Liao C.
        • Laporte A.D.
        • Spiegelman D.
        • Akçimen F.
        • Joober R.
        • Dion P.A.
        • Rouleau G.A.
        Transcriptome-wide association study of attention deficit hyperactivity disorder identifies 6 associated genes and phenotypes.
        Nat Commun. 2019; 10: 4450
        • Huckins L.M.
        • Dobbyn A.
        • Ruderfer D.M.
        • Hoffman G.
        • Wang W.
        • Pardiñas A.F.
        • et al.
        Gene expression imputation across multiple brain regions provides insights into 9 schizophrenia risk.
        Nat Genet. 2019; 51: 659-674
        • Hyde C.L.
        • Nagle M.W.
        • Tian C.
        • Chen X.
        • Paciga S.A.
        • Wendland J.R.
        • et al.
        Identification of 15 genetic loci associated with risk of major depression in individuals of European descent.
        Nat Genet. 2016; 48: 1031-1036
        • Hashimoto T.
        • Maekawa S.
        • Miyata S.
        IgLON cell adhesion molecules regulate synaptogenesis in hippocampal neurons.
        Cell Biochem Funct. 2009; 27: 496-498
        • Pischedda F.
        • Szczurkowska J.
        • Cirnaru M.D.
        • Giesert F.
        • Vezzoli E.
        • Ueffing M.
        • et al.
        A cell surface biotinylation assay to reveal membrane-associated neuronal cues: Negr1 regulates dendritic arborization.
        Mol Cell Proteomics. 2014; 13: 733-748
        • Sanz R.
        • Ferraro G.B.
        • Fournier A.E.
        IgLON cell adhesion molecules are shed from the cell surface of cortical neurons to promote neuronal growth.
        J Biol Chem. 2015; 290: 4330-4342
        • Milaneschi Y.
        • Simmons W.K.
        • Rossum EFC van
        • Penninx B.W.
        Depression and obesity: Evidence of shared biological mechanisms.
        Mol Psychiatry. 2019; 24: 18-33
        • Ochedalski T.
        • Subburaju S.
        • Wynn P.C.
        • Aguilera G.
        Interaction between oestrogen and oxytocin on hypothalamic-pituitary-adrenal axis activity.
        J Neuroendocrinol. 2007; 19: 189-197
        • Puder J.J.
        • Freda P.U.
        • Goland R.S.
        • Wardlaw S.L.
        Estrogen modulates the hypothalamic-pituitary-adrenal and inflammatory cytokine responses to endotoxin in women.
        J Clin Endocrinol Metab. 2001; 86: 2403-2408
        • Soares C.N.
        • Zitek B.
        Reproductive hormone sensitivity and risk for depression across the female life cycle: A continuum of vulnerability?.
        J Psychiatry Neurosci. 2008; 33: 331-343
        • Behl C.
        Oestrogen as a neuroprotective hormone.
        Nat Rev Neurosci. 2002; 3: 433
        • Hernández I.
        • Rosende-Roca M.
        • Alegret M.
        • Mauleón A.
        • Espinosa A.
        • Vargas L.
        • et al.
        Association of TMEM106B rs1990622 marker and frontotemporal dementia: Evidence for a recessive effect and meta-analysis.
        J Alzheimers Dis. 2015; 43: 325-334
        • Nicholson A.M.
        • Rademakers R.
        What we know about TMEM106B in neurodegeneration.
        Acta Neuropathol (Berl). 2016; 132: 639-651
        • Rutherford N.J.
        • Carrasquillo M.M.
        • Li M.
        • Bisceglio G.
        • Menke J.
        • Josephs K.A.
        • et al.
        TMEM106B risk variant is implicated in the pathologic presentation of Alzheimer disease.
        Neurology. 2012; 79: 717-718
        • Van Deerlin V.M.
        • Sleiman P.M.A.
        • Martinez-Lage M.
        • Chen-Plotkin A.
        • Wang L.-S.
        • Graff-Radford N.R.
        • et al.
        Common variants at 7p21 are associated with frontotemporal lobar degeneration with TDP-43 inclusions.
        Nat Genet. 2010; 42: 234-239
        • Byers A.L.
        • Yaffe K.
        Depression and risk of developing dementia.
        Nat Rev Neurol. 2011; 7: 323-331
        • Li J.Q.
        • Tan L.
        • Wang H.F.
        • Tan M.S.
        • Tan L.
        • Xu W.
        • et al.
        Risk factors for predicting progression from mild cognitive impairment to Alzheimer’s disease: A systematic review and meta-analysis of cohort studies.
        J Neurol Neurosurg Psychiatry. 2016; 87: 476-484
        • Ownby R.L.
        • Crocco E.
        • Acevedo A.
        • John V.
        • Loewenstein D.
        Depression and risk for Alzheimer disease: Systematic review, meta-analysis, and metaregression analysis.
        Arch Gen Psychiatry. 2006; 63: 530-538
        • Ichikawa T.
        • Baba H.
        • Maeshima H.
        • Shimano T.
        • Inoue M.
        • Ishiguro M.
        • et al.
        Serum levels of TDP-43 in late-life patients with depressive episode.
        J Affect Disord. 2019; 250: 284-288
        • Wainberg M.
        • Sinnott-Armstrong N.
        • Mancuso N.
        • Barbeira A.N.
        • Knowles D.A.
        • Golan D.
        • et al.
        Opportunities and challenges for transcriptome-wide association studies.
        Nat Genet. 2019; 51: 592
        • Finucane H.K.
        • Reshef Y.A.
        • Anttila V.
        • Slowikowski K.
        • Gusev A.
        • Byrnes A.
        • et al.
        Heritability enrichment of specifically expressed genes identifies disease-relevant tissues and cell types.
        Nat Genet. 2018; 50: 621

      Linked Article

      • Dissecting Genetically Regulated Gene Expression in Major Depression
        Biological PsychiatryVol. 89Issue 6
        • Preview
          One of the major aims of complex trait and population genetics is to identify causal genes underlying disease onset and progression. Causal genes may then be used to identify disease mechanisms and subsequently inform the development of more effective drugs and treatment strategies. Large-scale genome-wide association studies (GWASs), which test millions of single nucleotide polymorphisms (SNPs) for an association with a disease outcome, have greatly expanded the discovery of genetic risk loci underlying complex diseases.
        • Full-Text
        • PDF