|Home | About | Journals | Submit | Contact Us | Français|
Major depressive disorder (MDD) affects around 350 million people worldwide; however, the underlying genetic basis remains largely unknown. In this study, we took into account that MDD is a gene-environment disorder, in which stress is a critical component, and used whole-genome screening of functional variants to investigate the ‘missing heritability’ in MDD. Genome-wide association studies (GWAS) using single- and multi-locus linear mixed-effect models were performed in a Los Angeles Mexican-American cohort (196 controls, 203 MDD) and in a replication European-ancestry cohort (499 controls, 473 MDD). Our analyses took into consideration the stress levels in the control populations. The Mexican-American controls, comprised primarily of recent immigrants, had high levels of stress due to acculturation issues and the European-ancestry controls with high stress levels were given higher weights in our analysis. We identified 44 common and rare functional variants associated with mild to moderate MDD in the Mexican-American cohort (genome-wide false discovery rate, FDR, <0.05), and their pathway analysis revealed that the three top overrepresented Gene Ontology (GO) processes were innate immune response, glutamate receptor signaling and detection of chemical stimulus in smell sensory perception. Rare variant analysis replicated the association of the PHF21B gene in the ethnically unrelated European-ancestry cohort. The TRPM2 gene, previously implicated in mood disorders, may also be considered replicated by our analyses. Whole-genome sequencing analyses of a subset of the cohorts revealed that European-ancestry individuals have a significantly reduced (50%) number of single nucleotide variants compared with Mexican-American individuals, and for this reason the role of rare variants may vary across populations. PHF21b variants contribute significantly to differences in the levels of expression of this gene in several brain areas, including the hippocampus. Furthermore, using an animal model of stress, we found that Phf21b hippocampal gene expression is significantly decreased in animals resilient to chronic restraint stress when compared with non-chronically stressed animals. Together, our results reveal that including stress level data enables the identification of novel rare functional variants associated with MDD.
Major depressive disorder (MDD) causes considerable morbidity and mortality,1, 2, 3, 4, 5 and is a leading contributor to the global burden of disease.6 However, we know little about its underlying fundamental biology and the genes conferring susceptibility to this disorder.7
Decades of investigation have revealed little of the genetic basis of MDD. Last year, the first two loci for MDD were identified in Chinese women with severe symptoms, one near the SIRT1 gene and the other in an intron of the LHPP gene. However, neither gene has been replicated in European-ancestry populations,8 and a mega-analysis of several genome-wide association studies (GWAS) in MDD did not find any single nucleotide polymorphism (SNP) with genome-wide significance in European-ancestry populations.9 These results highlight the challenges facing this field despite the substantial inheritance of depression (37–38%),10, 11 which is even higher (~70%) when considering diagnostic unreliability.12 There have been various explanations postulated for this genetic conundrum, including disease heterogeneity and the types of genetic variation studied.13 Recently, the risk of depressive symptoms was associated with a rare missense variant in the LIPG gene; in addition, data from over 75000 and 230000 individuals with self-reported clinical diagnosis of depression and no self-reported history of depression, respectively, identified 15 genetic loci in European-ancestry populations.14, 15
In MDD, as in other mental disorders, there is a need for new genomic approaches to identify the ‘missing heritability’, which has arisen from applying the ‘common disease, common variant’ hypothesis. An alternative approach proposes that genetic risk can be better ascertained by considering nonlinear interactions between rare genetic variants that have stronger phenotypic effects.16, 17, 18, 19, 20 Coding and non-coding ‘functional’ SNPs (that is, those resulting in amino acid, splicing, regulatory or epigenetic changes) may have a major impact on phenotype.21 Indeed, important fractions of the current missing heritability of complex traits could be due to genetic interactions.20 Following this lead, we previously analyzed a small set of nonsynonymous SNPs and reported that their interactions with clinical and environment variables were significantly associated with MDD.22 Those findings provided a rationale to focus on functional SNPs.
MDD is a gene-environment disorder;3 however, genetic studies have not accounted for the fact that stressful events are a critical factor when selecting a control cohort. Unique features of this study include: (1) a cohort of cases and controls comprised of Mexican-Americans of the greater Los Angeles area who were mainly recent immigrants born in Mexico23 and who experienced significant levels of hyperactivation of the hypothalamic–pituitary–adrenal axis related to challenges, distress and acculturation issues related to immigration;24, 25 (2) type of genetic variations studied; (3) an ethnically diverse replication cohort with stress level data in controls. We performed whole-genome screening specifically of functional variants to investigate whether our experimental design could contribute to explain the ‘missing heritability’ of MDD. We obtained whole-genome sequencing for a small subset of individuals from both cohorts, and also performed an animal study to investigate the functional relevance of our genetic finding.
We used a Los Angeles Mexican-American cohort of 399 subjects aged 19–65 years: 203 with mild to moderate MDD cases (50.88%) and 196 controls (49.12%). There were no differences in age between MDD cases and controls (Supplementary Table S1). Participants gave written informed consent, and detailed demographic, epidemiological and clinical descriptions of this sample have been presented elsewhere.22, 23, 26 Briefly, our Mexican-American participants had at least three grandparents born in Mexico. They participated in a pharmacogenetic study with antidepressant drugs, and were assessed by the Structured Clinical Interview (SCID) for the Diagnostic and Statistical Manual of Mental Disorders (DSM), using the DSM-IV diagnosis of current, unipolar major depressive episode and a HAM-D21 (21-Item Hamilton Depression Rating Scale) score of 18 or greater with item number 1 (depressed mood) rated 2 or greater. This study was approved by the Institutional Review Boards of the Universities of California Los Angeles and Miami, USA, and by the Human Research Ethics Committees of the Australian National University and Bellbery, Australia, and it was registered in ClinicalTrials.gov (NCT00265291).22, 23, 26
Mexican-American MDD patients had comprehensive psychiatric and medical assessments in their primary language using diagnostic and rating instruments fully validated in English and Spanish. Exclusion criteria included active medical illnesses that could be etiologically related to the ongoing depressive symptoms, active suicidal intent, pregnancy/lactation, illicit drug use and/or alcohol abuse in the last 3 months, other major mental disorders, except for anxiety related disorders.22, 23, 26 Control age- and gender-matched Mexican-American individuals were recruited from the same Mexican-American community in Los Angeles and were in general good health but were not screened for medical or psychiatric illnesses.22, 23, 26
We used a European-ancestry cohort of 972 individuals (473 (51%) cases and 499 (49%) controls) for our replication study. The effects of age differences between MDD cases and controls were controlled in the linear mixed-effect models (LMEMs) maximization by introducing age as a covariate (Supplementary Tables S1 and S2). Participants gave written informed consent and were recruited under two protocols: (1) Mood disorder studies Münster (comprised of the Münster neuroimaging and the moodinflame studies, which have been conducted by the Department of Psychiatry and Psychotherapy, University of Münster, Münster, Germany), and (2) the Cognitive function and mood disorders study (conducted by the Discipline of Psychiatry, University of Adelaide, South Australia, Australia). The SCID/MINI was used to ascertain that healthy controls were free from lifetime history of psychiatric disorders; for this cohort, the main diagnostic and mood assessment instruments used were also DSM-IV criteria for MDD and HAM-D21, respectively. Previous traumatic events were assessed using the Childhood Trauma Questionnaire/Family Inventory of Life Events in cases and controls. Exclusion criteria included any neurologic abnormalities, substance-related disorders, psychotic/mania or hypomania symptoms, treatment with benzodiazepine, previous electroconvulsive therapy, and usual magnetic resonance imaging (MRI)-contraindications. European-ancestry cohort samples were studied under approved Human Research Ethics Committees protocols at the University of Münster, Germany, and University of Adelaide and Flinders University, South Australia, Australia.
Power estimation was done using the pwr27 package in R.28 Let n1 and n2 be the number of cases and controls, p1 and p2 the allele frequency for cases and controls respectively, and h the Cohen’s effect size, where . From this expression, it follows that and . Provided h and feasible values for p1 and p2, the difference in allele frequency between cases and controls can be calculated. Cohen29 suggests that h values of 0.2, 0.5 and 0.8 represent small, medium and large effect sizes, respectively.
Our analysis suggested that 200 cases and 200 controls would be sufficient to detect 80% true positives and a medium size effect (defined by the Cohen’s h parameter; Supplementary Figure S1) when m=100000 variants are tested for association (a value that overcomes the final number of variants used in the association analysis). On the basis of the effect sizes for the m=19 variants associated with MDD in the Mexican-American sample (Table 1A), the post-hoc power ranges between >60% (exm2249659, h=0.335) and >99% (exm1508600, h=0.643). Overall, these figures indicate that our study had power to detect medium to large effect sizes regardless of the strategy used (that is, GWAS or sequential-based).
Both cohorts were genotyped by the Australian Genome Research Facility (www.agrf.org.au) using the Illumina HumanExome BeadChip-12v1_A, which exonic content consists of >250000 markers representing diverse populations and a range of common conditions. Samples with calls below the Illumina (Illumina Australia and New Zealand, Scoresby, VIC, Australia) expected 99% SNP call rates were excluded. An individual was duplicated to test genotyping reliability and quality. The identity by descent (IBD) matrix between all pairs of individuals was estimated after linkage disequilibrium (LD) pruning and used for quality control and for the mixed linear models analyses.
GenomeStudio data were analyzed by SNP & Variation Suite (SVS) 7.6.7, Golden Helix’s (Golden Helix, Bozeman, MT, USA (http://www.goldenhelix.com). Parameters for excluding markers from analyses included: (i) deviations from Hardy–Weinberg equilibrium with P-values <0.05/m (where m is the number of markers included for analysis); (ii) a minimum genotype call rate of 90% (iii) the presence of more than two alleles; and (iv) monoallelism. The minor allele frequency (MAF) of 0.01 was the criterion for defining common (0.01) and rare variants (<0.01).
We filtered against functional prediction information available in the dbNSFP_NS_Functional_Predictions annotation track GRCh_37 to define variation with potential functional effect.30 This filter uses SIFT, PolyPhen-2, MutationTaster, Gerp++ and PhyloP.31, 32, 33 We used the SVS 7.6.7, Golden Helix’s Variant Classification module to examine the interactions between variants and gene transcripts to categorize variants based on their potential functional effects. Non-coding variants were tested for being harbored in splicing sequences, enhancers, CTCF (CCCTC-binding factor) binding sites, transcription factor binding sites, open chromatin regions, CpG islands, DNA and histone methylation sites, polymerase binding, and target sites of microRNAs and RNA-binding proteins.
We studied the association of MDD to functional variations using single- and multi-locus LMEMs34 with up to 10 steps in the backward/forward optimization algorithm. These models include both fixed (genotype markers, sex and years of education) and random effects (family or population structure), the latter to account for potential inbreeding by including the IBD matrix (which was estimated between all pairs of individuals after LD pruning in our analysis).34 A single-locus LMEM assumes that all loci have a small effect on the trait, whereas a multi-locus LMEM assume that several loci have a large effect on the trait.34 Both types of models were implemented in SVS 8.3.0. (Golden Helix). The optimal model was selected using a comprehensive exploration of multiple criteria including the Extended Bayes Information Criteria (eBIC), the Modified Bayes Information Criteria (mBIC) and the Multiple Posterior Probability of Association (mPPA). After the estimation process using the forward/backward algorithm was finished, the coefficients were extracted and a hypothesis test of the form H0,i: βi=0 vs H1,i: βi≠0 was performed for the ith common exonic functional variants to obtain the corresponding P-value (i=1,2,…,m). Thus, the collection P1, P2,…,Pm of P-values were corrected for multiple testing using the false discovery rate (FDR)35 and a method based on extreme-values theory.36 Because hypotheses testing were of the same type, correction was only performed on the resulting m P-values.35, 36 For the single-locus models, we estimated the Genomic Control ‘inflation factor’ λ, to evaluate potential stratification effects.
We used the regression- and the permutation-based kernel-based adaptive cluster (KBAC) methods to analyse rare exonic functional variants.18, 37 KBAC, implemented in the Golden Helix’s SVS 8.3.0, catalogs rare variant data within each of a number of regions/transcripts into multi-marker genotypes and determines their association with the phenotype, weighting each multi-marker genotype by how often that genotype was expected to occur according to control and MDD cases data and the null hypothesis that there is no association between that genotype and the case/control status.18, 37 Thus, genotypes with high sample risks are given higher weights that potentially separate causal from non-causal genotypes. A one-sided test was applied due to the weighting procedure and the P-values were estimated using 10000 permutations.
As an instrument of cross-validation we applied a sequential strategy of analysis by examining 50% of the Mexican-American individuals randomly chosen, and using the remaining set of individuals to replicate only those markers reaching genome-wide FDR <0.05 during the first round.
We used a replication European-ancestry cohort; functional variations harbored in genes, or nearby 30 kb, significantly associated with MDD in the Mexican-American sample were filtered in and tested for GWAS analysis of common and rare variations. For rare variants we performed KBAC with permutation testing on CCDS (consensus coding sequence) genes 15, UCSC (rfs://data.goldenhelix.com:80/rfs/CCDSGenes15-UCSC_2014-02-16_GRCh_37_Homo_sapiens.tsf:1). A matrix of IBD for the LMEMs analyses was generated with common variation filtered out after the targeting of genes associated in the first tier. The criterion for replication was a P-value of 0.10 after correction for multiple comparisons using FDR.
Haplotype analyses were performed using Golden Helix. In the Mexican-American and European-ancestry cohorts, we interrogated the chromosomal region with the strongest replication signals to characterize LD blocks and transmission. Lewontin’s D statistics was applied to those analyses because it is unaffected by rare allele frequencies.38
A small subset of both cohorts (15 Mexican-American samples and 10 Australians of European-ancestry samples) was sequenced using Illumina HiSeq 2000 (BGI, Shenzhen, Guangdong, China) or HiSeq X (Garvan Institute, Sydney, NSW, Australia). Whole-genome sequencing paired-end reads were aligned to the human reference genome (hg19, Genome Reference Consortium GRCh37) using Burrows-Wheeler Aligner (BWA)39 to get SAM (sequence alignment/map) format files. Then SAMtools40 was used to convert the SAM to the BAM (Binary version of a SAM file) format files. After internal sorting, BAM files were merged into one BAM file (where one sample may include several lanes of sequencing reads). We used the mpileup command in SAMtools to collect summary information from the BAM file, calculate the likelihood of data given each possible genotype, and store the likelihoods in a binary file. We then piped the output to SAMtools/BCFtools,41 which makes the SNV (single nucleotide variant)/INDEL (small insertions and deletions) calling to generate the VCF (variant call format) files. After that, we used ANNOVAR42 to annotate the SNP/INDEL information and their classification details in those VCF files.
We evaluated potential common ontogenetic and cellular function processes of the genes disclosed by the GWAS analysis on rare and common variants using Metacore 6.8 software build 29806 (GeneGo, St. Joseph, MI, USA).
We used the UK Brain Expression Consortium (UKBEC) web based server Braineac (Brain eQTL Almanac, http://www.braineac.org/) to understand whether some genetic polymorphisms have significant statistical association with transcription level (expression quantitative trait loci, eQTL) in ten brain regions (cerebellum, frontal cortex, hippocampus, medulla, occipital cortex, putamen, substantia nigra, thalamus, temporal cortex and central white matter).43
Procedures were approved by the Animal Ethics Committees of the South Australian Health and Medical Research and Flinders University and are in accordance with the Australian Code for the Care and Use of Animals for Scientific Purposes (8th edition, 2013). Virus- and antibody-free young adult male Sprague–Dawley rats (150–200g) obtained from Charles River (Margate, Kent, UK) were individually housed in Green Line IVC Sealsafe PLUS cages (Tecniplast, Varese, Italy) in a temperature- (22±1°C) and light- (12h cycles, lights on at 07:00 hours) controlled, stress-free and specific pathogen-free environment with water and standard regular chow ad libitum. They were allowed to habituate at least 5 days before the initiation of experimental procedures. After the baseline behavioral testing, rats were assigned to one of two groups so that the range of floating time was similarly distributed and there was no floating time differences between groups before treatment: (i) chronic restraint stress (CRS, n=27) and (ii) non-CRS (n=10) groups. Moreover, the CRS group was further classified into CRS resilient and CRS non-resilient subgroups (see ‘Statistical analyses’ section below). The investigator was not blinded to group allocation; however, most of the behavioral data were collected by a camera coupled to a software (EthoVision, Noldus Information Technology, Wageningen, The Netherlands). The number of animals was decided based on pilot studies carried out in the lab.
Flat-bottom clear acrylic restraint containers (20.3 × 8.3 cm; cat no. 544-RR Plas Labs, Lansing, MI, USA) were used as previously described.44 For 14 consecutive days during daytime (09:00–16:00 hours) CRS rats were submitted to six consecutive hours of restraint, after which they were unrestrained and returned to their home cages. Non-CRS animals were submitted to behavioral testing but not restrained.
We used the automatic video-tracking EthoVision XT video tracking software (EthoVision) to record and analyze the behavior and activity of our animals during the forced swim test (FST). Tests were performed between 09:00 and 12:00 hours after a 60 min habituation to the testing room. At baseline and after the CRS paradigm, the pre-test (training) was carried out for 10 min and was followed 24 h later by the test. Rats were individually tested in a glass cylinder (45cm height and 30cm diameter), which contained 30 cm of water (the rat’s hindlimbs did not reach the cylinder’s floor) at 23°C for 5min. Activity was recorded by one perpendicular camera located in the front; automated measurements of floating/immobility (<12% of distance moved) and struggling/highly mobile activity (>18.5% of distance moved) were obtained twice per second using the EthoVision XT software.45 After the test, rats were dried and placed in a 30 °C drying environment for ~10–15 min. FST was performed at baseline and four weeks later, at which time the CRS paradigm had been completed.
After the last behavioral testing, animals were allowed to rest at least 2 days before being killed. To avoid the confounding effects of circadian rhythms, animals were killed between 10:00 and 12:00 hours. Brains were remove from the skull, immediately submerged in RNAlater solution and stored overnight at 4 °C; supernatant was then removed and brains were dissected and stored at −80°C until use. Total RNA was isolated from hippocampi using the Purelink RNA Mini Kit (Life Technologies, Mulgrave, VIC, Australia), and DNase digestion was performed using the Purelink DNase set (Life Technologies). RNA was quantified using a spectrophotometer (NanoDrop 2000, Thermo-Fisher Scientific, Waltham, MA, USA) and 500 ng total RNA was reversed transcribed to cDNA using iScript RT supermix (BioRad, Hercules, CA, USA) and random hexamer primers.
Primers for each gene of interest were exon spanning, designed using the IDT primer quest tool (Integrated DNA Technologies, Baulkham Hills, NSW, Australia) or the primer premier 5 software (Premier Biosoft International, Palo Alto, CA, USA). Amplicon lengths were between 76 and 135 bp. We employed the Bestkeeper46 and the geNorm47 softwares to select 2 out of 5 reference genes we tested. The following primer pairs were used:Phf21bF: 5′-CAGCGGAAGGCCTTAAAGAA-3′ Phf21bR: 5′-CACTGTCTTGTGGGTGACATAG-3′ Rps18F: 5′-TTCAGCACATCCTGCGAGTA-3′ Rps18R: 5′-TTGGTGAGGTCAATGTCTGC-3′ GapdhF: 5′-CCATTCTTCCACCTTTGATGCT-3′ GapdhR: 5′-TGTCATACCAGGAAATGAGCTTCA-3′.
A standard curve of pooled, serially diluted cDNA was run for the Phf21b gene and for housekeeping genes (Gapdh and Rps18) using the QuantStudio 7 Flex Real-Time PCR system (Thermo-Fisher Scientific). cDNA samples were diluted 1:4 and run in triplicates. Primer sets were tested for optimal dissociation curves with amplification efficiencies between 94 and 107%. To check for genomic contamination, a minus-reverse transcriptase control of each sample was run in an quantitative real-time reverse transcriptase (RT-qPCR) experiment. The geometric mean of both housekeeping genes was used to calculate the results. Only reactions with threshold cycle (CT) standard deviation values 0.3 were accepted.
Mirroring our clinical approach, our analyses took into consideration the stress level displayed by animals under non-chronically stressed conditions, as our rats were individually housed and social isolation can increase anxiety- and depressive-like behaviors.48, 49, 50 In both groups, we excluded non-chronically stressed animals that displayed increased floating time at baseline (in the CRS group) or increased averaged floating time (in the non-CRS group), which was defined as floating time higher than the post-CRS FST floating time mean. We classified the CRS group into CRS resilient (below average baseline and post-CRS floating times) and CRS non-resilient. Differences between groups were analyzed by the Student-paired t-test or one-way ANOVA when appropriate. The significance level for each of these effects was set at P<0.05. We used the Pfaffl method to calculate RT-qPCR data.51 Statistical analyses were performed using GraphPad Prism 6 (GraphPad Software, La Jolla, CA, USA).
We obtained Illumina HumanExome BeadChip-12v1_A genotype data from 399 Mexican-American subjects (203 MDD cases and 196 controls) and 972 European-ancestry subjects (499 MDD cases and 473 controls). After filtering out markers not meeting either quality control criteria or variability requirements 83898 variants remained for analyses (Figure 1). In the Mexican-American cohort, a total of 19 common SNPs in 18 genes were significantly associated with MDD at the genome-wide FDR <0.05 (Figure 2; Table 1A); two SNPs were found in the MUC5B gene. Sixteen nonsynonymous SNPs out of the nineteen genome-wide associated SNPs predicted functional protein effects and eleven MDD-associated SNPs were rare variants in controls (MAF<0.01). Haplotype analysis reported the clustering of markers with significant association to four conspicuous regions harboring the OR2T12 (Chr 1), TMRM150B (Chr 19), SLX4 (Chr 16) and TRPM2 (Chr 21) genes (Supplementary Figures S2A–D).
We next applied the KBAC method to 47296 rare exonic variants with potential functional effects. Twenty-seven genes were significantly associated with MDD in the Mexican-American cohort (Table 1B). Functional SNPs harbored in five of these genes, ANO8, CNTNAP1, EMR2, HOMER3, UNC13D, were significantly associated with MDD in common and rare variants analyses (Tables 1A and and1B1B).
For cross-validation, we performed a sequential strategy by first analyzing data of 50% randomly chosen Mexican-American individuals, and then using the remainder 50% to replicate only those markers reaching genome-wide FDR <0.05 during the first round. Ten of the thirteen SNPs significantly associated with MDD in the sequential analyses (Table 1C) also achieved significance (FDR<0.05) in the GWAS analysis (Table 1A).
Functional variations harbored in or 30kb around genes significantly associated with MDD in the Mexican-American sample were next tested by GWAS analysis of common and rare variations in the European-ancestry cohort. The rare variant analysis using the KBAC method identified the PHF21B gene (Table 2), which is adjacent to the PRR5-ARHGAP8 gene within the chromosome 22q.13 region identified in Mexican-Americans. It is noteworthy to mention that the frequencies of genetic variants were different between both cohorts (see ‘WGS analyses’ section below).
Haplotype analyses were performed in Chr 22q13.1 in both cohorts encompassing the region showing the highest replication signals in Table 2, namely the PRR5, ARHGAP8 and PHF21B genes (Supplementary Figure S3). The LD structure of this chromosomal region appears to be more conserved in the European-ancestry cohort than in the Mexican-American cohort. It is evident that this chromosomal region is transmitted as a heterogeneous block. Therefore, diversity of the results in our two cohorts might be explained by the complex evolutionary heterogeneity of this region.
We performed WGS to better understand variations within and between each of the two cohorts we studied. WGS analyses showed that Mexican-American individuals had significantly more SNVs (50% more) and less INDELs (small insertions and deletions) when compared with Australian individuals of European-ancestry (Table 3; Supplementary Table S3). Our data are compatible with findings from the International HapMap 3 Consortium and the 1000 Genomes Project Consortium, as they respectively showed that individuals with African ancestry carry an increased number of variants and rare variants, and within Europeans the Spanish population carry excess of rare variants.52, 53 As the HapMap Mexican-American samples from Los Angeles and our Mexican-American cohort were recruited from the same location by our bilingual research team, their median ancestry proportions are estimated to be 45% Indigenous American, 49% European and 5% African.54
MDD-associated genes in the Mexican-American cohort were overrepresented in the following Gene Ontology (GO) processes: (1) natural immune response (granuloma formation, natural killer cell degranulation, synaptic vesicle priming, natural killer cell activation, germinal center formation), (2) G-protein-coupled glutamate receptor signaling pathway, (3) detection of chemical stimulus involved in smell sensory perception, (4) cellular response to stimulus, (5) regulation of macromolecule metabolic processes, (6) muscle contraction, (7) cell surface receptor process (G-protein coupled receptor signaling pathway, system process, neurological system process), (8) regulation of ion transmembrane transport (regulation of calcium ion transport and activity) and (9) positive regulation of cellular metabolic process (Supplementary Table S4).
Data extracted from the Braineac web server (UK Brain Expression Consortium, UKBEC) revealed that PHF21B gene variants significantly change eQTL in all ten listed brain areas (Supplementary Table S5). Seven listed variants change hippocampus eQTL (Table 4).
Our genetic analyses were performed on two cohorts that, uniquely, took into consideration levels of stress data in control groups to minimize the impact of individuals genetically at-risk for MDD in the control population and thereby maximize the likelihood of identifying novel MDD-associated variants. This approach revealed that rare functional variants in the PHF21B gene are associated with MDD in both cohorts and brain eQTL showed PHF21B gene variants could significantly change expression in several brain areas. To further analyze the association between PHF21B, depression and stress, we exposed male Sprague–Dawley rats to CRS (see ‘Materials and methods’ section), and analyzed the effect on floating/immobility time in the FST, which is a measure of behavior despair and depressive-like behavior. CRS resulted in small but significantly increased floating time in the FST in CRS non-resilient animals but not in the CRS resilient ones (Figure 3a; Supplementary Table S6). Importantly, Phf21b mRNA levels were significantly decreased in the CRS-resilient group when compared with the non-CRS group, but not in the CRS non-resilient group (Figure 3b). Therefore, Phf21b gene expression level may reflect a state of resilience to chronic stress.
We presented results that support the assumption that MDD is a syndrome of considerable genetic heterogeneity, as we identified many common and rare functional variants that confer susceptibility to MDD in the Mexican-American cohort. Those variations withstood the use of stringent correction measures related to genetic stratification and repeated testing. Interestingly, 11 out of 19 SNPs with significant genome-wide association to MDD consisted of rare variants in the control population; the data strongly support a Mendelian type of inheritance for many of the rare variants. Remarkably, five genes were identified in both common and rare variant analyses, and one gene had two SNPs that were significantly associated with MDD in our GWAS analysis.
Our replication strategy in the European-ancestry cohort included testing functional variations harbored in genes, and extended to include nearby 30 kb, significantly associated with MDD, in GWAS analysis of common and rare variations, in the Mexican-American. This strategy replicated the PHF21B (PHD finger protein 21B, also known as BHC80L or PHF4) gene, which encodes the 531 amino acid PHD finger protein 21B. The PHF21B gene and the neighboring ARHGAP8 (Rho GTPase activating protein 8) gene provided the strongest signals in our replication analysis; however, the ARHGAP8 gene did reach the significance threshold (Table 2). Furthermore, as the ARHGAP8 and PRR5 (Proline rich 5) genes were initially thought to be one gene because of readthrough transcripts,55 this chromosomal locus requires further scrutiny.
The PHF21B gene may be a putative tumor suppressor gene whose loss of function results from reduced expression by hypermethylation or gene loss.56 This gene is expressed in the brain and its variants affect its expression in several brain regions, including the frontal cortex and the hippocampus (Brain eQTL Almanac, http://www.braineac.org/); though, its central nervous system functions are unknown. However, its prominent paralog, the PHF21A gene, encodes BHC80 (a component of a BRAF35/histone deacetylase complex), which is a transcriptional repressor during neurodevelopment highly expressed in the brain.57, 58, 59 Both the PHF21A and the PHF21B genes have the PHD zinc finger domain, the DNA-binding domain and the region required for transcriptional repression, and localize in the nucleus.56 The PHF21B gene is located in Chr 22q13.31, which is a genomic region associated with the Phelan-McDermid syndrome (or 22q13.3 deletion syndrome/22q13.3DS), a rare neurodevelopmental disorder that typically presents with generalized developmental delay, intellectual disability, delayed speech and seizures, and involves the SHANK3 gene.60 The SHANK3 gene encodes a protein that has a role in synapse formation and dendritic spine maturation, and mutations in this gene cause autism spectrum disorder and schizophrenia.61, 62 It is noteworthy that the nearby 22q11.2 deletion syndrome [22qDS/velocardiofacial syndrome (VCFS)/DiGeorgio syndrome/CATCH22] is one of the most common multiple anomaly syndromes in humans and it also includes among its manifestations higher rates of psychiatric disorders, including mood and schizophrenia.63, 64 A large percentage (79%) of 22qDS children and adolescents have one psychiatric diagnosis and 12.5–14% meet DSM criteria for MDD.65, 66, 67
Results of our haplotype and sequential analyses support that the TRPM2 gene was significantly associated with MDD; thus, these data may be considered as providing a replication for the association of MDD and the TRPM2 gene, which was previously associated by positional candidate approach with susceptibility to bipolar disorder and unipolar disorder.68 The transient receptor potential cation channel, subfamily M, member 2, TRPM2, is a member of the melastatin-related transient receptor channel family that is involved in oxidative stress-induced cell death and inflammation processes and has a role in endotoxin-provoked cytokine production; TRPM2-mediated calcium influx influences the reactive oxygen species-induced signaling cascade responsible for chemokine production and exacerbates inflammation via the NLRP3 (NLR family, pyrin domain containing 3) inflammasome.69
Pathway analysis offered insight into the complexity of mood regulation and MDD stereotypical behaviors and symptoms. Many of the uncovered enriched processes are novel or underexplored, such as cellular response to stimulus, muscle system process, regulation of metabolic process, sensory perception of chemical stimulus, regulation of ion transmembrane transport and protein–DNA complex disassembly, and others have been active areas of investigation as complementary/alternative hypothesis for the classical monoamine hypothesis of depression, such as neuroimmune mediation/inflammation and glutamate receptor signaling pathway and growth/positive regulation of macromolecule biosynthetic process.70, 71, 72, 73, 74, 75 The olfactory system and its central connections may have a role in affective behavior, and bilateral olfactory bulbectomy has been used as a rodent model for depression since the early 1980’s.76, 77 Our data may help narrow the focus of future investigations toward specific aspects of these broad investigative areas.
Our animal study provided evidence that hippocampal Phf21b gene expression modulates the response to chronic stress as CRS-resilient animals had decreased hippocampal Phf21b mRNA levels. However, the lack of non-stressed group housed animals may have restricted our ability to understand the baseline level of this gene in non-socially isolated animals. GO enrichment analysis indicated that MDD-associated genes are involved in regulation of metabolic processes; therefore, food and water restriction during the CRS paradigm could have impacted post-stress gene expression levels. To minimize the effects of food and water restriction in our experiments, CRS was performed during the light phase when rodents do not usually consume a significant amount of food or water, as they are nocturnal animals.
MDD is clearly a gene-environment disorder,3 but most genetic studies have not accounted for stressful life events in the control population. These critical environmental factors have been used in the present work to help minimize the inclusion of genetically susceptible individuals in the control sample. Our Mexican-American cohort is comprised of first-generation individuals (60%)23 who have experienced significant levels of stress and hyperactivation of the hypothalamic–pituitary–adrenal axis related to acculturation issues.24, 25 In contrast to the Mexican-American controls represented a group of individuals who were highly resilient to significant levels of enduring stress, in our replication cohort the stress levels in controls were also included in our analyses. Furthermore, as most of the variations reported here in the Mexican-American cohort were rare, each gene contained low-frequency mutations at many different sites, which could help explain the challenges experienced in this field, as the power of genome-wide association studies in MDD have been greatly reduced because the critical assumption that there were little allelic heterogeneity within loci may not apply;78 this assumption may have been harder to negate in populations of European-ancestry due to their drastically reduced number of SNVs (Table 3; Supplementary Table S3). Our findings suggest that the ‘missing heritability’ in MDD could be at least partly explained by rare variants.
In summary, we identified common and rare variations in a total of 44 genes that may confer susceptibility to MDD in a Mexican-American cohort. Most of these variations were rare and resulted in amino acid (that is, likely functional) changes. Replication of the PHF21B gene in an ethnic diverse population and the finding that Phf21b gene expression modulates the chronic stress response in rats corroborate the strength of our findings. Our findings also provide a set of common and rare genetic variants associated with MDD that require replication. Our clinical cohorts were small; therefore, further replication studies are warranted with larger cohorts recruited under a similar study design to replicate and/or identify additional genes associated with MDD. As PHF21B gene expression could be modulated by methylation,56 future studies should investigate whether the methylation status in this locus is influenced by MDD diagnosis or stress, as epigenetics may account for part of the missing heritability in complex traits.79 No previous studies have looked at Phf21b levels in stress/depressive-like behavior; therefore, additional animal studies are also needed to systematically characterize hippocampal Phf21b gene expression level in non-stressed groups at baseline, and during various types of acute and chronic stress paradigms.
We have been supported by grants APP1051931 and APP1070935 (MLW), and APP1060524 (BB) from the National Health and Medical Research Council (Australia), the German Research Foundation (UD, Grant FOR 2107, DA1151/5-1), NIH Grant GM61394 (JL and MLW), the German Australian Institute for Translational Medicine (SRB and JL), and institutional funds from the Australian National University and the South Australian Health and Medical Research Institute. We are grateful for the contributions of Israel Alvarado, Deborah L. Flores, Rita Jepson, Lorraine Garcia-Teague, Patricia Reyes, Isabel Rodriguez, Gabriela Marquez and the University of California, Los Angeles General Clinical Research Center (UCLA GCRC) staff in recruiting and caring of the Mexican-American participants.
Supplementary Information accompanies the paper on the Molecular Psychiatry website (http://www.nature.com/mp)
The authors declare no competing interest and no income from pharmaceutical companies; an intellectual property application has been prepared to include the genetics findings of this work.