Introduction
Recurrent major depressive disorder (MDD) is a complex phenotype with limited established associations with biological markers and molecular genetic pathways. The China, Oxford and Virginia Commonwealth University Experimental Research on Genetic Epidemiology (CONVERGE) study and the Psychiatric Genomics Consortium (PGC) have made advances in the molecular genetic etiology of MDD, confirming that MDD is highly polygenic, with aggregation of top loci effects accounting for a small proportion of the variance in case–control status (Peterson et al., Reference Peterson, Cai, Bigdeli, Li, Reimers, Nikulova, Webb, Bacanu, Riley, Flint and Kendler2017; Major Depressive Disorder Working Group of the PGC et al., Reference Wray, Ripke, Mattheisen, Trzaskowski, Byrne, Abdellaoui, Adams, Agerbo, Air, Andlauer, Bacanu, Baekvad-Hansen, Beekman, Bigdeli, Binder, Blackwood, Bryois, Buttenschon, Bybjerg-Grauholm, Cai, Castelao, Christensen, Clarke, Coleman, Colodro-Conde, Couvy-Duchesne, Craddock, Crawford, Crowley, Dashti, Davies, Deary, Degenhardt, Derks, Direk, Dolan, Dunn, Eley, Eriksson, Escott-Price, Kiadeh, Finucane, Forstner, Frank, Gaspar, Gill, Giusti-Rodriguez, Goes, Gordon, Grove, Hall, Hannon, Hansen, Hansen, Herms, Hickie, Hoffmann, Homuth, Horn, Hottenga, Hougaard, Hu, Hyde, Ising, Jansen, Jin, Jorgenson, Knowles, Kohane, Kraft, Kretzschmar, Krogh, Kutalik, Lane, Li, Li, Lind, Liu, Lu, MacIntyre, MacKinnon, Maier, Maier, Marchini, Mbarek, McGrath, McGuffin, Medland, Mehta, Middeldorp, Mihailov, Milaneschi, Milani, Mill, Mondimore, Montgomery, Mostafavi, Mullins, Nauck, Ng, Nivard, Nyholt, O'Reilly, Oskarsson, Owen, Painter, Pedersen, Pedersen, Peterson, Pettersson, Peyrot, Pistis, Posthuma, Purcell, Quiroz, Qvist, Rice, Riley, Rivera, Saeed Mirza, Saxena, Schoevers, Schulte, Shen, Shi, Shyn, Sigurdsson, Sinnamon, Smit, Smith, Stefansson, Steinberg, Stockmeier, Streit, Strohmaier, Tansey, Teismann, Teumer, Thompson, Thomson, Thorgeirsson, Tian, Traylor, Treutlein, Trubetskoy, Uitterlinden, Umbricht, Van der Auwera, van Hemert, Viktorin, Visscher, Wang, Webb, Weinsheimer, Wellmann, Willemsen, Witt, Wu, Xi, Yang, Zhang, Arolt, Baune, Berger, Boomsma, Cichon, Dannlowski, de Geus, DePaulo, Domenici, Domschke, Esko, Grabe, Hamilton, Hayward, Heath, Hinds, Kendler, Kloiber, Lewis, Li, Lucae, Madden, Magnusson, Martin, McIntosh, Metspalu, Mors, Mortensen, Muller-Myhsok, Nordentoft, Nothen, O'Donovan, Paciga, Pedersen, Penninx, Perlis, Porteous, Potash, Preisig, Rietschel, Schaefer, Schulze, Smoller, Stefansson, Tiemeier, Uher, Volzke, Weissman, Werge, Winslow, Lewis, Levinson, Breen, Borglum and Sullivan2018). Thus, supporting examination of variants across the entire genome, rather than just top loci, will likely be most informative to etiology.
Statistical genomics can be used to examine both global polygenic risk and risk within established gene pathways. With genomic data, it is increasingly feasible to begin to isolate and examine molecular drug targets for MDD. One way to work toward this aim is by inspecting molecular pathway (MP) enrichment in relation to severe MDD diagnosis, or in relation to symptom presentation in severe cases. It is also important to present caveats: MP analyses are fraught with potential for spurious findings and require (1) large sample sizes, (2) strict correction for (a) multiple comparison and (b) average background polygenic signal, (3) a grain of salt, given the diffuse and quite polygenic signal from even very large GWAS meta-analyses, and thus (4) replication. It is also important to note that available database information about MPs is far from complete.
Keeping these important points in mind, we also know that many MPs are enriched for MDD signal, and that medication targeting-specific neurotransmitter pathways are efficacious. Thus, it is possible that MPs have some influence on the presentation of and risk for severe MDD, and potentially on drug response. Given the aim of genomic research to identify promising drug targets for MDD, the identification of potential molecular classifiers of drug response is highly desirable.
To this aim, we examined global polygenic risk for MDD, and also polygenic risk for MDD specific to MPs, in the presentation of severe MDD in a large case–control sample of Han Chinese descent (N = 10 060) (CONVERGE Consortium, 2015). CONVERGE aimed at identifying genetic risk factors for recurrent MDD among a rigorously ascertained cohort, and presents several advantages for genetic study. The sample is entirely female, thus reducing heterogeneity from sex differences. Further, unlike most studies of its size, CONVERGE subjects all have four Han Chinese grandparents, thus reducing noise stemming from population stratification. Finally, to increase severity and potential genetic signal, all subjects were ascertained for recurrent depression, rather than a single episode. This ascertainment method, and the size of the study, allowed CONVERGE to identify and replicate risk loci for MDD (CONVERGE Consortium, 2015).
This analysis used summary statistics from a leave-one-batch-out GWAS of MDD within CONVERGE. We constructed Bayesian subject-level genome-wide MDD polygenic risk scores (PRS) and MDD PRS specific to pruned MPs. Scores and covariates were regressed onto case status, and pathway scores with significant signal above and beyond global MDD polygenic risk were then examined by within-pathway, gene-based PRS.
Methods and materials
Sampling
Recurrent MDD cases and healthy controls were recruited from 51 mental health centers and psychiatric departments of general medical hospitals across 21 provinces of China. Details have been previously published (CONVERGE Consortium, 2015). We controlled for potential clinical heterogeneity by recruiting only female participants, and to control for ethnic stratification, only participants whose grandparents (all four) were of Han Chinese descent were recruited to participate. Cases and controls [age M (s.d.) = 44.4 (8.9) and 47.7 (5.6), respectively] were excluded for diagnosis of bipolar disorder, any psychotic symptoms occurring outside of the context of a major depressive episode, and any diagnosis of intellectual disability. Description of the clinical features of the cases can be found in Edwards et al. (Reference Edwards, Docherty, Moscati, Bigdeli, Peterson, Webb, Bacanu, Hettema, Flint and Kendler2018). Cases had at least two major depressive episodes with the first episode occurring prior to age 50, and could not have abused drugs or alcohol prior to their first episode of depression. Controls were clinically screened to rule out prior depressive episodes and had to be at least 40 years of age, past the age window of typical MDD onset. The study protocol was approved centrally by the Ethical Review Board of Oxford University, and by the ethics committees of all of the participating hospitals in China. All participants provided written informed consent.
DNA sequencing
DNA extraction, sequencing, genotyping, and imputation details have been reported (Cai et al., Reference Cai, Bigdeli, Kretzschmar, Li, Liang, Hu, Peterson, Bacanu, Webb, Riley, Li, Marchini, Mott, Kendler and Flint2017). Briefly, DNA was extracted from saliva using Oragene and sequenced reads were obtained from Illumina Hiseq machines aligned to Genome Reference Consortium Human Build 37 patch release 5 (GRCh37.p5) with Stampy (v1.0.17) with default parameters (Lunter and Goodson, Reference Lunter and Goodson2011). Reads consisting of base quality ⩽5 or containing adaptor sequencing were filtered out. Alignments were indexed in BAM format (Lunter and Goodson, Reference Lunter and Goodson2011) using Samtools (v0.1.18) (Li et al., Reference Li, Handsaker, Wysoker, Fennell, Ruan, Homer, Marth, Abecasis and Durbin2009), and PCR duplicates were marked for downstream filtering using Picardtools (v1.62). The Genome Analysis Toolkit's (GATK, version 2.6) BaseRecalibrator (McKenna et al., Reference McKenna, Hanna, Banks, Sivachenko, Cibulskis, Kernytsky, Garimella, Altshuler, Gabriel, Daly and DePristo2010) created recalibration tables to screen known SNPs and INDELs in the BAM files from dbSNP (version 137, excluding all sites added after version 129). GATKlite (v2.2.15) was used for subsequent base quality recalibration and removal of read pairs with improperly aligned segments as determined by Stampy.
Calling and imputation of genotypes
GATK's UnifiedGenotyper VariantRecalibrator (version 2.7-2-g6bda569) was used on post-BSQR files for variant discovery and genotyping at all polymorphic SNPs in 1000 G Phase1 ASN panel (The 1000 Genomes Project Consortium et al., Reference Abecasis, Auton, Brooks, DePristo, Durbin, Handsaker, Kang, Marth and McVean2012), as well as variant quality score recalibration. A sensitivity threshold of 90% was applied for imputation after optimizing for transition to transversion ratios. Genotype likelihoods were calculated using a binomial mixture model implemented in SNPtools (version 1.0) (Wang et al., Reference Wang, Lu, Yu, Gibbs and Yu2013), and imputation was performed at sites with no reference panel using BEAGLE (version 3.3.2) (Browning and Browning, Reference Browning and Browning2007). A second round of imputation was performed at biallelic polymorphic SNPs using 1000G Phase 1 ASN haplotypes as a reference panel. To determine the final number of SNPs, we applied a conservative inclusion threshold for SNPs: (1) a p value for violation HWE >10–6, (2) information score >0.9, and (3) minor allele frequency in CONVERGE >0.5%.
Diagnostic assessments
Participant interview sessions lasted approximately 2 h. Interviewers were largely trained psychiatrists with a small number representing post-graduate medical students or psychiatric nurses, and all were clinically trained by the CONVERGE team for at least one week. Interviews were recorded and included an assessment of psychopathology, demographic characteristics, and psychosocial functioning. Trained editors listened to a portion of the interviews to provide ratings of interview quality. An assessment of lifetime adversity (binary yes/no) reflected a composite of endorsements on self-reported stressful life event and childhood sexual abuse measures (Peterson et al., Reference Peterson, Cai, Dahl, Bigdeli, Edwards, Webb, Bacanu, Zaitlen, Flint and Kendler2018). In CONVERGE, adversity has been observed to contribute significant variance in genetic signal, and was included as a feature in our follow-up analyses only to determine whether genetic signal is influenced, in any direction, by presence of lifetime adversity. We excluded participants who had incomplete assessment information or were lacking high-quality genetic data, to arrive at a final 4728 controls and 5612 case samples for analysis.
Polygenic risk for MDD and pathway-based risk as a predictor of case status
MDD PRSs were constructed globally (genome-wide) and within each of the MPs, using summary statistics from leave-one-batch-out GWAS analyses of the CONVERGE sample. A list of pathways was derived from independent, previously published pathway enrichment statistics from the PGC (Schizophrenia Working Group of the Psychiatric Genomics Consortium, 2014). These pathways were deemed to likely be the most relevant to psychiatric disorders broadly. Gene lists corresponding to each pathway were extracted from the relevant gene ontology database, with the number of genes per pathway ranging from 8 to 598 (mean = 9). Python-based LDpred (Vilhjalmsson et al., Reference Vilhjalmsson, Yang, Finucane, Gusev, Lindstrom, Ripke, Genovese, Loh, Bhatia, Do, Hayeck, Won, Kathiresan, Pato, Pato, Tamimi, Stahl, Zaitlen, Pasaniuc, Belbin, Kenny, Schierup, De Jager, Patsopoulos, McCarroll, Daly, Purcell, Chasman, Neale, Goddard, Visscher, Kraft, Patterson and Price2015), a polygenic risk scoring software, allows for the modeling of LD to weight the relative contributions of syntenic variants to the outcome phenotype. LDpred uses postulated proportions of causal variants in the genome as Bayesian prior probabilities for PRS calculations, and a range of eight different priors were used (proportions of 1, 0.3, 0.1, 0.03, 0.01, 0.003, and 0.001) to construct scores. All pathway scores were LD pruned and merged based on excessive correlation between scores (R 2 < 0.5). Stepwise regressions were run with all PRS variables and critical covariates as determined by the CONVERGE Consortium (2015). MDD was treated as a dichotomous variable (despite likely heterogeneity) in order to avoid predicting additional clinical information within cases such as symptom count.
We tested two models for the relationship between PRSs and MDD status. Model 1 included all risk scores and the two primary ancestry principal components. Model 2 included all risk scores, primary ancestry principal components, and the presence/absence of lifetime adversity (a significant moderator of the genetic signal in previous CONVERGE GWAS). Follow-up analyses included similar associations of any significant pathway-based PRS associations with age at onset (AAO), family history of MDD, and number of lifetime episodes within the MDD cases.
Results
For the primary analysis of this paper, i.e. finding pathways associated with MDD status, we employed a stepwise regression that pared down the number of pathways PRSs to only the significant ones (after adjusting for multiple testing). Statistics for the best-fitting model are presented in Table 1, each parameter in the model remaining significant after FDR correction. The best-fitting stepwise regression included three robust predictors: Global MDD, one ancestry PC, and the MDD PRS for the GO:0017144 drug metabolism pathway. Model 2 results are presented in Table 2. Global MDD risk and GO:0017144 pathway-based risk remained significant in Model 2, where lifetime adversity was entered as a covariate (and was found to be the most robust predictor). Of note, the stepwise models captured no increase in pathway-based risk by genome-wide MDD PRS. Standardized effect sizes, presented as squared Pearson partial correlation coefficients in both Tables 1 and 2, were obtained after adjusting for multiple testing (Bigdeli et al., Reference Bigdeli, Lee, Webb, Riley, Vladimirov, Fanous, Kendler and Bacanu2016). We mention that these analyses appear to be well powered, given that (even after adjusting for multiple testing) we estimated to have a power >80% to detect PRS explaining only 0.003 of variation in MDD liability, which is half of the in-sample estimates for GO:0017144.
Best-fitting stepwise regression model predicting major depressive disorder case status from global major depressive disorder polygenic risk score, molecular pathway-based polygenic risk score, and ancestry
MDD, major depressive disorder; L1O, leave-one-out cross-validation scoring procedure; PRS, polygenic risk score; GO, gene ontology database.
Best-fitting stepwise regression model accounting for adversity
MDD, major depressive disorder; L1O, leave-one-out cross-validation scoring procedure; PRS, polygenic risk score; GO, gene ontology database.
Lifetime adversity is binary (1 = lifetime adversity and 0 = no lifetime adversity).
While not the main goal of this work, due to their likely reduced power, we use a gene PRS-type secondary/exploratory analysis to suggest genes in the drug metabolism pathway that might be most important to MDD. Of the 33 genes included in the pathway, most of which reflect protein coding for cytochrome P450 and flavin containing monooxygenase (https://icgcportal.genomics.cn/genesets/GO:0017144), 88 803 variants on nine genes passed variant filtering and LD pruning. Exploratory gene-based analysis included MDD PRS calculation specific to the nine genes involved in the GO:0017144 pathway. Exploratory LDpred analyses of the PRSs associated with the 33 genes yielded three genes with suggestive signals (Table 3). CYP2C19 and CBR1 both yielded FDR q-values of 0.05. Follow-up analyses of the GO:0017144 pathway in the MDD cases resulted in a significant association with AAO in the expected direction, after correction for multiple testing (t = 3.388, p = 0.0007 without global MDD PRS in the model; t = 3.382, p = 0.0008 with global MDD PRS in the model). Gene-level analyses of CBR1 and CYP2C19 predicting AAO were suggestive but not significant (p < 0.08).
Suggestive genes in drug metabolism pathway after false discovery rate correction
GO, gene ontology database; Chr:Mbp, chromosome and position; FDR, false discovery rate corrected p value.
Discussion
In the CONVERGE sample, MDD polygenic risk specific to the GO:0017144 pathway significantly predicted MDD case status above and beyond genome-wide polygenic risk for MDD. Exploratory analyses suggest that among genes in this pathway, CYP2C19 and CBR1 appear to contribute the most to MDD liability, as their effect on trait remained significant after accounting for both MDD polygenic risk and for an established genetic risk moderator, lifetime adversity. These scores were based on MDD GWAS weights entirely within-ancestry, and did not rely on summary statistics from Northern European ancestry GWAS. This is the first study to test MDD PRS specific to MPs in relation to recurrent MDD, and preliminary results suggest that GO:0017144 may be relevant to recurrent MDD and to clinical presentation.
MDD is genetically complex, and it is important to try to identify areas of focus on the genome with which to inform drug target research. These results are potentially relevant to pharmacogenomics. Several common antidepressants (tricyclic antidepressants and SSRIs, including citalopram, escitalopram, sertraline, amitriptyline, and clomipramine) are substrates for CYP2C19 (Sim et al., Reference Sim, Nordin, Andersson, Virding, Olsson, Pedersen and Ingelman-Sundberg2010). More importantly, CYP2C19 possesses epoxygenase activity and metabolizes various prostaglandins (arachidonic acid, linoleic acid, eicosapentaenoic acid) to active epoxide products, which have vasodilatory, anti-inflammatory and antinociceptive effects (Fabbri et al., Reference Fabbri, Tansey, Perlis, Hauser, Henigsberg, Maier, Mors, Placentino, Rietschel, Souery, Breen, Curtis, Lee, Newhouse, Patel, O'Donovan, Lewis, Jenkins, Weinshilboum, Farmer, Aitchison, Craig, McGuffin, Schruers, Biernacka, Uher and Lewis2018). Similarly, CBR1 plays a protective role in oxidative stress, neurodegeneration, and apoptosis. CBR1 normally breaks down prostaglandin E2, which may mediate the aversive/negative affect component of inflammatory pain (Singh et al., Reference Singh, Zajdel, Mirrasekhian, Almoosawi, Frisch, Klawonn, Jaarola, Fritz and Engblom2017).
Because MDD is likely quite heterogeneous, it is important to try to identify genetic classes of cases that differ phenotypically (Milaneschi et al., Reference Milaneschi, Lamers, Peyrot, Baune, Breen, Dehghan, Forstner, Grabe, Homuth, Kan, Lewis, Mullins, Nauck, Pistis, Preisig, Rivera, Rietschel, Streit, Strohmaier, Teumer, Van der Auwera, Wray, Boomsma and Penninx2017). This view is supported by research attempting to genetically stratify severe psychopathology (Howard et al., Reference Howard, Clarke, Adams, Hafferty, Wigmore, Zeng, Hall, Gibson, Boutin, Hayward, Thomson, Porteous, Smith, Murray, Haley, Deary, Whalley and McIntosh2017). Recent research by Howard et al. has successfully replicated a genetic classification of MDD using PRS of ~20 medical and psychiatric phenotypes. Future research would benefit from examination of whether high-risk PRS classes are phenotypically distinct (whether clinical features differ) and also whether they are distinct with respect to risk within drug-relevant MPs.
Were there enough non-European GWAS to complete a similar stratified analysis to CONVERGE, we would have attempted this in conjunction with our analysis of MP-based genetic MDD risk. PRS for several medical and psychiatric phenotypes can currently be calculated for European samples based on European GWAS, but unfortunately the field lacks sufficient East Asian GWAS for analysis in CONVERGE. We hope that future GWAS of non-European samples will facilitate additional meaningful insights.
Strengths
A strength of this study is that ascertainment successfully minimized genetic heterogeneity due to sex differences and to ancestry. Moreover, ascertainment was limited to recurrent depression rather than a single episode, potentially increasing the strength of the genetic signal. In addition, this large cohort provided sufficient power to detect even small effects, and deep phenotyping allowed for the additional inspection of lifetime adversity in relation to genetic signal. We included a limited number of pathways – only those implicated in severe psychopathology – in order to limit Type I error.
Limitations
Strict ascertainment of females with four Han Chinese grandparents may limit the generalizability of GO:0017144 pathway findings. Because the study rigorously ascertained recurrent MDD, and not single episode, it is also possible that the effects observed here are limited to severe and/or recurrent MDD. One source of potential bias stems from the inclusion of the pathways implicated in severe psychopathology in Europeans rather than in individuals of Han Chinese ancestry. In addition, we did not detect signal for some candidate pathways (e.g. serotonergic or noradrenergic) that we might have expected to be relevant to MDD. Null results could stem from insufficient coverage across enough genes in the pathway after quality control, or pathway signal being obscured by the global polygenic signal for MDD that we accounted for in our models. This would be consistent with the small genetic effect sizes we find in such a large MDD cohort.
Finally, there are problems inherent to any study of gene pathways. These include the field's very limited molecular understanding of gene pathways and related sequelae, incompleteness of pathway databases, conflicting opinions about the independence/interdependence of gene pathways, and the elevated polygenicity of all psychiatric traits. To the first three points, we assert that it is studies like this that further develop the pathway literature and inform pathway database efforts. To the fourth point, it is critical for studies to account for genome-wide polygenic risk when testing for any signal relevant to a small portion of the genome. These points underscore the importance of independent replication.
Author ORCIDs
Anna R. Docherty, 0000-0001-7139-7007.
Acknowledgements
The authors gratefully acknowledge the support of all partners in hospitals across China, without whom this research would not be possible.
Financial support
This work was supported by the Wellcome Trust (J.F., grant numbers WT090532/Z/09/Z, WT083573/Z/07/Z, WT089269/Z/09/Z); the National Institutes of Health (A.D., grant number K01MH093731; A.E., grant number AA021399; D.A., grant number KMH109765; K.K., grant number MH100549; R.P., grant number K01MH113848; S.B., grant numbers MH100560, AA022717); the Simons Foundation (A.D., grant number 513631); the Brain & Behavior Research Foundation (A.D.); and the American Foundation for Suicide Prevention (A.D.).
Conflict of interest
None.
