INTRODUCTION

Major depression is the world's most common mental illness and a major source of disability (Ustun et al, 2004). Depending on its longitudinal course, major depression can appear as single or recurrent episodes (unipolar major depressive disorder, MDD) or alternating with episodes of mania or hypomania (bipolar disorder, BD). The latest estimate of the lifetime prevalence of mood disorders (MD) in Western Europe is around 14% (Alonso et al, 2004), and similar figures have been reported in the United States of America (Kessler et al, 2003).

MD patients commonly show biological rhythm-related symptoms, such as characteristic disturbances in the sleep/wake cycle, diurnal mood changes, and a periodic pattern of symptom recurrence and remission. In addition, alterations in the circadian pattern of core body temperature and neuroendocrine secretion have been documented (reviewed by McClung, 2007a). Furthermore, the so-called ‘clock’ genes, which have a central role in the generation and control of circadian rhythms, are also involved in the regulation of the monoaminergic pathways (Hampp et al, 2008; McClung, 2007b), whose implication in the pathophysiology and treatment of MD is well known (Belmaker and Agam, 2008; Nestler and Carlezon, 2006). In fact, the antidepressant fluoxetine alters the expression of clock genes in mouse brain regions that have been implicated in mood regulation, such as the hippocampus and striatum (Uz et al, 2005). There is also evidence that fluoxetine exerts its neuronal effects by influencing circadian timekeeping of serotonergic tone (Sprouse et al, 2006). Moreover, lithium, a mood-stabilizing drug used as a first-line treatment of BD and also as an augmentation agent of antidepressant treatment in MDD, is known to modulate circadian rhythms (Abe et al, 2000). What is more, non-pharmacological therapeutic interventions based on the regulation of biological rhythms, such as sleep deprivation, light therapy, and interpersonal and social rhythm therapy, have proved effective in the treatment of MD (Frank et al, 2000; Terman and Terman, 2005; Wu and Bunney, 1990).

In mammals, circadian rhythms in physiology and behavior are generated and synchronized by the endogenous cellular clock located in the suprachiasmatic nucleus of the hypothalamus. This master clock is regulated at the molecular level by complex mechanisms involving positive and negative transcriptional/translational feedback loops, which drive the circadian rhythmicity of clock gene transcripts (reviewed by Dardente and Cermakian, 2007; see Supplementary Figure S1 for a detailed description).

Although circadian abnormalities have long been implicated in the pathophysiology of MD (Healy, 1987), and although it has been hypothesized that genetic variations in clock-related genes are involved in the liability to MD (Bunney and Bunney, 2000; Mitterauer, 2000), the functional and molecular elements on which these effects are based have not been conclusively determined to date. Some reports suggested evidence for association of genetic variation in ARNTL, PER3, TIMELESS, and NR1D1 genes with BD (Kishi et al, 2008; Mansour et al, 2006; Nievergelt et al, 2006; Severino et al, 2009), and a significant interaction between three circadian genes (BHLHB2, CSNK1E, and CLOCK) has also been associated with BD (Shi et al, 2008). Far less is known regarding the involvement of circadian genes in the vulnerability to MDD. No evidence of association was found between a CLOCK gene variation in MDD (Bailer et al, 2005; Desan et al, 2000; Kishi et al, 2009a), although a positive association with fluvoxamine response was reported in Japanese MDD patients (Kishi et al, 2009b).

There is some evidence implicating the genetic variation in circadian genes with clinical and course variables in MD patients. For instance, the CLOCK gene has been associated with age at onset (Benedetti et al, 2004), illness recurrence (Benedetti et al, 2003), evolution of insomnia (Serretti et al, 2005), actimetric sleep and diurnal activity patterns (Benedetti et al, 2007), and information processing (Benedetti et al, 2008b). Moreover, the PER3 gene has been linked with MD features such as age at onset, response to SSRI treatment, circadian mood oscillations, and characteristics of temperament (Artioli et al, 2007; Benedetti et al, 2008a). Finally, the PER2, ARNTL, and NPAS2 genes have been associated with seasonal affective disorder (Johansson et al, 2003; Partonen et al, 2007).

All these findings suggest that genes controlling circadian rhythms are likely to be involved in susceptibility to MD (either unipolar MDD or BD). To explore this possibility, we ran a population-based genetic case–control study and genotyped tagSNPs covering the genomic region of 19 candidate circadian genes, some of which have not been previously investigated in MD: ADCYAP1, ARNTL, ARNTL2, BHLHB2, BHLHB3, CLOCK, CRY1, CRY2, CSNK1E, DBP, NPAS2, NR1D1, PER1, PER2, PER3, RORA, TIMELESS, VIP, and VIPR2.

MATERIALS AND METHODS

Subjects and Clinical Assessment

The clinical sample consisted of 534 unrelated Caucasoid patients of Spanish origin with MD (335 MDD and 199 BD) and 440 unrelated controls. Patients were consecutively recruited from two tertiary centers in Barcelona, Bellvitge University Hospital (HUB, n=445) and Hospital de la Santa Creu i Sant Pau (HSP, n=89) between 2004 and 2007 (female: 67.4%; mean age at recruitment: 55.4 years, SD 15.5; mean age at onset of illness: 38.5 years, SD 15.3). The patients were diagnosed by two experienced psychiatrists through a structured interview (First et al, 1997) according to the Diagnostic and Statistical Manual of Mental Disorders IV criteria (APA, 1994) for Mood Disorders. All patients had a history of at least one major depressive episode of moderate severity. Depression severity was assessed using the 21 item HAM-D administered by a psychiatrist during the depressive index episode (Hamilton, 1960). Exclusion criteria were age under 18 years, additional past or present psychiatric diagnosis other than MD, past or present history of psychoactive substance abuse except for nicotine, and severe medical disease. Patients were eligible for the study when both research examiners agreed on all criteria. In case of doubt or disagreement, a third psychiatrist was consulted and a consensus decision was reached. Clinical assessment included information on sociodemographic and clinical variables. Other variables such as age at onset and previous episodes were obtained by means of a personal interview, supported by the available medical information.

The control group consisted of 440 unrelated participants randomly selected from population registers, who were interviewed at Primary-care Centers in the same geographical area as the referral hospitals. Their mean age was 43.1 years (range: 18–75). Screening for psychiatric diseases was performed according to the 28-scaled Global Health Questionnaire (GHQ) (Lobo et al, 1986) and an interview conducted by a physician to provide a full medical history and details on lifestyle and habits. Diseases were recorded according to the International Classification of Diseases Version 9 (ICD9). Only unrelated citizens with no personal history of psychiatric disease, no clinically documented first-degree family history of psychiatric disease, and a GHQ score <7 were considered for inclusion in the study. The percentage of excluded participants was gender-independent. The ethical review boards of the corresponding recruitment centers approved the study. All subjects, both cases and controls, gave written informed consent to participate.

Marker Selection and Genotyping

We selected 248 single-nucleotide polymorphisms (SNPs) at 19 circadian genes (Supplementary Table S1 and Supplementary Figure S1): adenylate cyclase-activating polypeptide 1 (ADCYAP1); aryl hydrocarbon receptor nuclear translocator-like (ARNTL); aryl hydrocarbon receptor nuclear translocator-like 2 (ARNTL2); basic helix-loop-helix domain containing class B, 2 and 3 (BHLHB2 and BHLHB3); circadian locomoter output cycles kaput protein (CLOCK); cryptochrome 1 and 2 (CRY1 and CRY2); casein kinase 1, epsilon (CSNK1E); D site of albumin promoter (albumin D-box)-binding protein (DBP); neuronal PAS domain protein 2 (NPAS2); nuclear receptor subfamily 1, group D, member 1 (NR1D1); period 1, 2, and 3 (PER1, PER2 and PER3); RAR-related orphan receptor A (RORA); timeless (TIMELESS); vasoactive intestinal peptide (VIP) and vasoactive intestinal peptide receptor 2 (VIPR2). To select tagSNPs from the genomic region where each of the candidate genes are located, we first downloaded the publicly available genotypes corresponding to CEU trios of European descent of the HapMap project dataset (Phase II data freeze, assembly NCBI b35, dbSNP b125). Only SNPs with a unique mapping location and a minor allele frequency (MAF) higher than 5% were considered. Bins of common SNPs in strong LD, as defined by r2 higher than 0.85, were identified within this data set by the use of HapMap-LDSelect-Processor, which uses the ‘LD Select’ method (Carlson et al, 2004) to process HapMap genotype dump format data corresponding to the locus region. This approach led to the selection of 215 tagSNPs, but we also included 24 additional SNPs due to their previously reported involvement in MD or circadian phenotypes in the literature, and 9 non-synonymous SNPs which, although not present in HapMap, could be functionally interesting (Supplementary Tables S1 and S2). SNPs were genotyped using the SNPlex Genotyping System (Applied Biosystems, Foster City, CA), a multiplexed assay format that allows the simultaneous genotyping of up to 48 SNPs. We submitted the selected SNPs to the SNPlex design pipeline and selected the ones which, according to the design rules of the technique, were expected to develop good assays. The high-throughput genotyping assays were performed at the CeGen genotyping facilities, in the Barcelona Node (Centro Nacional de Genotipado, Genoma España) following the manufacturer's protocol for the SNPlex Genotyping System (Applied Biosystems) genotyping process. GeneMapper software v3.7 (ABI, Foster City, CA) was further used to analyze the raw data and automatically assign individual genotypes. All genotyping was performed blind with respect to phenotype. Details of the quality control of genotypes and genotyping for population admixture are given in Supplementary Methods section.

Statistical Analysis

Raw genotype data were available for 534 cases and 440 controls, resulting in a final data set of 974 individuals. Each polymorphism was tested in controls to ensure the fitting with Hardy–Weinberg equilibrium (HWE). We removed the SNPs with call rates of less than 90% (n=29), monomorphic in our population (n=6), or outside HWE in the control subset (n=5; p<0.01). Thus, 39 SNPs were removed (one SNPs was both out of HWE and had a call rate under 90%) and the final analyses were performed with 209 available SNPs (Supplementary Table S2). The genetic heterogeneity between centers was ruled out by comparing the allelic frequencies of the patients from the two hospitals (HUB and HSP) with the Haploview v.4.1 software (Barret et al, 2005), which carries out χ2-tests and allows permutation tests (set to 10 000) to correct for multiple testing.

We first carried out an exploratory association analysis in the total MD sample to achieve the maximum power for detecting common variants of modest effect, and we then performed the association analyses after splitting the total MD sample into unipolar MDD and BD subsamples, as previously performed (Kishi et al, 2009d). QuantoVersion 1.2.3 was used for power calculations (Gauderman and Morrison, 2006; Gauderman, 2002).

Allelic case–control association analysis was performed using χ2-tests as implemented in Haploview software v 4.1. Multivariate methods based on logistic regression analyses were used to perform the genotype association analyses. We estimated the odds ratios (OR) adjusted by sex and their 95% confidence intervals (95% CI). The SNP effects were fitted under four models of inheritance: codominant, dominant, recessive, and a log-additive model. The best model was selected using the Akaike information criteria. These analyses were performed using the SNPassoc R package (Gonzalez et al, 2007). Statistical significance was assessed by a permutation procedure to estimate the significance of the best result (1000 permutations) at a gene-wide level (a within-gene corrected p-value for each of the 19 genes tested); thereafter, these gene-wide p-values were corrected experiment-wide, also with a permutation approach, for all the SNPs tested across the 19 genes. The correction for the allelic tests and haplotype analysis was performed with the Haploview software forcing the permutation of individual SNPs along with the haplotypes in the examined blocks. As Haploview does not perform genotypic analyses, we permuted the corresponding p-values with the SNPassoc R package. Finally, considering that we carried out an exploratory analysis in the whole MD sample and then performed a stratified analysis on BD and MDD separately, permuted p-values below 0.025 were considered statistically significant.

The multiple-marker association analyses (also adjusted by sex) were restricted a priori to the SNPs that presented statistically significant associations with MD at gene-wide level, in order to minimize type I errors. All haplotype frequencies were estimated using the expectation maximization algorithm as implemented in the Haplo.stats package in the R programming language (Lake et al, 2003). In order to provide a comprehensive assessment of the haplotype subsets within the gene, we used a sliding window approach with the Haplo.stats R package. To evaluate potential additive effects between the SNPs that were statistically significant after correction for multiple comparisons identified in the single-marker association analysis at gene-wide level, regression nested models adjusted by sex were compared by a likelihood ratio test in the R programming language.

Comparison with Previous Whole-Genome Scan Data

Several groups have recently performed whole-genome association studies of BD (Baum et al, 2008; Ferreira et al, 2008; Sklar et al, 2008; WTCCC, 2007). Results for variants that showed positive association after correction for multiple comparisons at gene-wide level in our sample were only available from the WTCCC study, which genotyped 1900 BD type I patients from the United Kingdom and a shared set of 3000 controls (WTCCC, 2007). We therefore downloaded the data available online (http://www.wtccc.org.uk/info/summary_stats.shtml) corresponding to genotype counts and tested them for associations with the WTCCC's BD cases and their shared set of controls using the SNPassoc R Package (Gonzalez et al, 2007).

RESULTS

The demographic and clinical characteristics of the sample are summarized in Table 1. Power calculations determined that, under a log-additive model, our MD case–control sample had 80% power to detect a genotype-relative risk of 1.3 for a variant with a MAF of 24%.

Table 1 Demographic and Clinical Characteristics of the MD Sample (MDD+BD)
Full size table

We tested our cohort for population stratification calculating the genomic inflation factor based on the median χ2-test comparing the observed and expected distributions. We estimated that the inflation factor (λ) was equal to 1.00 (mean χ2-statistic is 1.006), its expected value being under no population stratification, thus indicating that overall population structure had negligible impact on the case–control association. Furthermore, a Q-Q plot (Supplementary Methods and Results) showed that observed p-values gradually depart from expected p-values at the extreme tail, a pattern that suggests that our association findings were more likely due to true genetic variation than other reasons such as genotyping errors, sample relatedness, or potential population stratifications. Finally, we found no evidence of genetic heterogeneity between the two groups from different hospitals in the same geographic area (HUB and HSP), as there were no significant differences in their allele frequencies after correcting by permutation (data not shown). All further analyses were performed using the combined data set.

Association Analyses and Evaluation of Additive Effects in MD

Single-marker analysis considering the combined MD data set (MDD plus BD) identified nominally statistically significant associations (uncorrected p<0.05) with one or more SNPs in 10 circadian genes: ADCYAP1, ARNTL, ARNTL2, CLOCK, CRY1, NPAS2, PER2, PER3, VIP, and VIPR2 (Supplementary Table S3). After permutation correction for the number of SNPs analyzed for each gene, variants located on the CRY1 (rs2287161), NPAS2 (rs11123857), and VIPR2 (rs885861) genes remained significantly associated with MD patients (Table 2). We also found a statistically significant additive effect of rs2287161 in CRY1, rs11123857 in NPAS2, and rs885861 in VIPR2: according to our results, carriers of the C allele in rs2287161, carriers of the G allele in rs11123857, and patients homozygous for the C allele in rs885861 were three times more likely to suffer from MD (OR=3.27; 95% CI: 2.1–5.1; p-value=0.005).

Table 2 Statistically Significant Results after Multiple Testing Correction* of the Single-Marker Association Study in Mood Disorders Adjusted by Sex
Full size table

The multimarker analyses revealed globally significant evidence for association in the case of CRY1 and NPAS2 genes. Particularly, the rs2287161–rs10861683 (CT) combination in CRY1 was significantly more frequent among MD patients (OR=1.80; 95% CI: 1.39–2.32; global permuted p-value=0.019), and the rs11123857–rs17662394 (GT) combination in NPAS2 was also over-represented among MD patients (OR=1.56; CI: 1.19–2.06; global permuted p-value=0.025) (Table 3).

Table 3 Statistically Significant Results of the Haplotype Analyses of CRY1 and NPAS2 Genes in Mood Disorder Patients and Controls
Full size table

Association Analyses and Evaluation of Additive Effects in MDD

Analyzing the MDD subsample separately, we found nominally statistically significant associations (p<0.05) with one or more SNPs in the following circadian genes: ADCYAP1, ARNTL, ARNTL2, CLOCK, CRY1, NPAS2, PER2, PER3, VIPR2 (also significant in the total MD sample) as well in DBP and NR1D1 (Supplementary Table S4). After correction for multiple testing at gene-wide level, only rs2287161 in CRY1 and rs11123857 in NPAS2 remained significantly associated with MDD (Table 2). After correction for multiple comparisons at experiment-wide level, only rs2287161 in CRY1 remained significant (p=0.007 under an additive model). Moreover, and in agreement with the results obtained with the total MD sample, a significant additive effect of the C allele in rs2287161 (CRY1) and the G allele in rs11123857 (NPAS2) was found in the MDD subsample (OR=2.38; 95% CI: 1.73–3.28; p-value=0.0009).

After applying a multiple comparison correction by permutation at gene-wide level, the multimarker analyses also showed evidence of association between the same haplotype in CRY1, being again over-represented among MDD patients (OR=1.91; CI: 1.42–2.57; p-value=0.013), although this association is not statistically significant at experiment-wide level (p=0.032) (Table 3).

Association Analyses and Evaluation of Additive Effects in BD

In the BD sample, we found nominally statistically significant associations (p<0.05) with one or more SNPs in 12 circadian genes: ADCYAP1, ARNTL, ARNTL2, CLOCK, CRY1, NPAS2, PER3, VIP, and VIPR2, as in the total MD sample, whereas BHLHB3, CSNK1E, and RORA also showed associations with BD (Supplementary Table S5). After correcting for the number of SNPs analyzed in each of the tested genes, only rs10462028 in CLOCK and rs17083008 in VIP remained significantly associated with BD (Table 2). Furthermore, we detected a significant protective additive effect against the presence of BD diagnosis when the subjects were heterozygous (AG) for rs10462028 (CLOCK) and also carried the G allele in rs17083008 (VIP) (OR=0.20; 95% CI: 0.09–0.43; p-value=0.003). The multimarker analysis did not reveal a haplotypic association within the CLOCK gene or the VIP gene in our BD subsample (data not shown).

Comparison with Previous Whole-Genome Scan Data in BD

We compared the results of our SNPs that were significant at the gene-wide level in any subgroup of our patients with those of the WTCCC's bipolar cohort. Specifically, we compared the following SNPs: rs10462028 in CLOCK, rs17083008 in VIP, rs2287161 in CRY1, rs11123857 in NPAS2, and rs885861 in VIPR2. First, we performed an association analysis in the WTCCC BD cases and their shared set of controls (58BC and UKBS), and observed an association with rs10462028 (CLOCK) under a log-additive model (OR=0.91; 95% CI: 0.84–1.00; crude p-value=0.04). When the two sets of WTCCC controls were considered separately, rs10462028 (CLOCK) was statistically significant for the UKBS set of controls under a dominant model (OR=0.85; 95% CI: 0.74–0.97; crude p-value=0.02). When we considered only the 58BC set of controls, we found a significant association with rs2287161 (CRY1) under the log-additive model (OR=0.89; 95% CI: 0.81–0.99; crude p-value=0.02), as in our samples. In spite of these nominal results, if we take into account that we analyzed five different SNPs, none of the tested variants would stand multiple testing correction.

DISCUSSION

In this study, we performed a case–control genetic association study of MD patients and 19 candidate genes involved in the regulation of circadian rhythmicity through the extensive genotyping of tagSNPs covering all the genomic regions where the genes are located. We first carried out an exploratory association analysis in the combined sample of major MD, as MDD and BD not only share features such as phenotypic clinical expression of major depressive episodes, abnormalities in circadian rhythms, pathophysiological pathways, or therapeutic issues, but there may also be an overlap in their genetic susceptibility (McGuffin and Katz, 1989; McGuffin et al, 2003). Moreover, a proportion of patients with MDD can develop manic or hypomanic episodes and become bipolar in the future (Akiskal et al, 1995) or even fit into de concept of ‘bipolar spectrum’ without showing obvious manic symptomatic polarity (Akiskal and Pinto, 2000). In addition, some clinical data support a unitary conceptualization of MD (Benazzi, 2007; Cassano et al, 2004). However, they differ in several clinical, pathophysiological, and therapeutic features (Ghaemi et al, 2004; Yatham et al, 1997) and independent loci have been reported to be associated with one or other diagnosis (Kato, 2007). This supports the idea that, although the disorders may have a similar genetic background, some genes may be specific to each one. To take this possibility into account and to facilitate comparisons with previous studies that considered only one of the diagnoses, we performed the association analyses considering MDD and BD as separate subsamples.

Association signals that remained significant after correction for multiple testing at the gene-wide level were observed for CRY1, NPAS2, and VIPR2 in the combined MD sample. In the stratified analysis based on polarity, the same SNPs and two-marker haplotype in CRY1 remained associated only in the MDD subgroup. In fact, our strongest association was an SNP located 3′ near CRY1 gene in MDD, which remained statistically significant after applying multiple testing correction at experiment level (p=0.00008). In spite of these apparently specific effects on the depressive pole in our sample, there is additional evidence for this association. First, CRY1 is located on chromosome 12q23–q24.1, a region that has been reported by several groups to be linked to BD (Curtis et al, 2003; Dawson et al, 1995; Degn et al, 2001; Detera-Wadleigh, 1999; Morissette et al, 1999); second, it was nominally associated in our BD sample (p=0.04 under a log-additive model and p=0.03 in the allelic association). Finally, the data from the WTCCC cohort of bipolar patients also revealed a nominal significant association with the same CRY1 SNP (p=0.02). It is reasonable to think that our weak statistical result in the BD subsample is probably due to a loss of statistical power after the sample was stratified, which reduced the size of the BD subset (although the results of the WTCCC study, with a much larger sample size, were not much stronger). Thus, the genetic effect of this gene on MD susceptibility is probably modest, and the correction for multiple testing that we applied led necessarily to the exclusion of potentially interesting results. Based on the above evidence, the involvement of CRY1 in BD cannot be ruled out.

The significant variant rs2287161 in the CRY1 gene is an intergenic SNP located 3′ downstream from the gene's polyadenylation site, which does not cause any obvious deleterious genetic change, thus suggesting the presence of 3′ downstream regulatory elements that are able to modulate CRY1 gene expression over long distances. On the other hand, rs11123857, located on intron 16 of NPAS2, is predicted to be an intronic enhancer and thus could cause an abnormal pre-mRNA splicing, which might modify expression levels. Finally, an interesting result of this study is that we detected the presence of genetic additive effects between these two SNPs. NPAS2, which is a homolog of CLOCK, forms transcriptionally active heterodimers with BMAL1, and CRY1 exerts its repression through the lock of the heterodimer in an active state (Dardente et al, 2007). So it is not surprising that variants in CRY1 or NPAS2 may each have a particular effect on the maintenance of molecular rhythmicity and that, when found together, they increase the effect of the same neurobiological network. On the other hand, changes in NPAS2 have been previously related with winter depression (Partonen et al, 2007) and may have a parallel effect in other cognitive and behavioral phenotypes. For instance, NPAS2 knockout mice have been reported to be deficient in complex emotional memory (Garcia et al, 2000), a trait whose neural substrate in humans is in the amygdala and hippocampal complex (Phelps, 2004), structures that have been implicated in the pathophysiology of MD (Blumberg et al, 2003; Sheline, 2003). Furthermore, emotional memory in humans is mainly enhanced across periods of sleep, and it has consequently been hypothesized that this sleep-dependent affective modulation could be related to the initiation and/or maintenance of mood disturbance in MD (Walker, 2009). Sleep disturbances are key features of depressive symptomatology, with subjective sleep complaints recorded in more than 80% of depressed patients (Reynolds and Kupfer, 1987). Therefore, variants in CRY1 and NPAS2 could have an impact on mood regulation through sleep mechanisms, as both genes seem to have an important role in sleep homeostasis (Franken et al, 2006; Wisor et al, 2002) and genetic variation in core clock genes is linked to circadian rhythm sleep disorders in humans (Ebisawa et al, 2001; Toh et al, 2001).

With respect to the analysis of the BD subset, CLOCK (rs10462028, in the 3′ downstream region of the gene) and VIP (rs17083008, in the promoter/regulatory region) were found to be disorder-specific. Evidence for association of an interaction between three circadian genes including CLOCK with BD has been previously reported (Shi et al, 2008), and, although not replicated, decreased levels of VIP in plasma and CSF have been found in lithium-treated euthymic bipolar patients (Berrettini et al, 1985). As only rs17083008 in VIP is predicted to have a potential functional role since it alters a transcription factor-binding site, it raises the possibility that alternative non-typed variants in LD with our significant SNPs are the true functional ones. In spite of this, it is interesting to note that our significant SNP in CLOCK (rs10462028) gave nominally significant results in the BD WGAs from the WTCCC study (p=0.02), supporting the association of CLOCK with BD in samples with a European genetic background.

The specific association of these genes with BD but not with MDD suggests their specific involvement in mania–hypomania phenotypes. Transgenic mice carrying a mutation in the Clock gene display an overall behavioral profile that is strikingly similar to human mania, including hyperactivity, decreased sleep, lowered depression-like behavior, and lowered anxiety. Moreover, chronic administration of the mood stabilizer lithium reverts many of these behavioral responses to wild-type levels (Roybal et al, 2007). Hyperactivity is also a key clinical feature of attention-deficit hyperactivity disorder (ADHD) and an association between a 3′-UTR SNP (rs1801260) and adult ADHD has been reported (Kissling et al, 2008). Interestingly, the same SNP was nominally associated in our BD sample (p=0.004 under a codominant model), suggesting the involvement of CLOCK gene variation in the susceptibility to clinical features shared by the mania pole of BD and ADHD, such as hyperactivity, inattention, and impulsivity. In recent years, there has been growing evidence of the critical involvement of VIP and its receptor VIPR2 in circadian timekeeping (Maywood et al, 2006; Vosko et al, 2007). Vip- and Vipr2-deficient mice exhibit disruptions in their ability to express circadian rhythm in constant conditions, together with a loss of synchrony in electrical activity and weakened light response (Colwell et al, 2003; Harmar et al, 2002). Therefore, it is believed that while 24-h transcriptional/translational feedback loops in individual neurons may ultimately control the expression of circadian rhythmicity, agents outside these loops, such as VIP, are also necessary (Loh et al, 2008). This would explain the genetic additive effect that we observed between CLOCK and VIP variants, as it is conceivable that disruption in both mechanisms can induce more profound effects in the resulting phenotype.

Regarding the overlap with previous positive findings from studies investigating the association of circadian genes in MD, and despite the limited power of our BD subsample, we found nominally statistically significant associations in ARNTL and PER3 genes, in agreement with previous results suggesting an association with BD in family based association designs (Mansour et al, 2006; Nievergelt et al, 2006). Two recent studies report nominal significant associations of genetic variation in NR1D1 in BD (Kishi et al, 2008; Severino et al, 2009) that are not replicated in our BD subset of patients, but conversely we found a positive nominal association in our MDD subsample, suggesting the possible involvement of NR1D1 in other MD-related phenotypes. However, NR1D1 gene does not seem to have a major role in the antidepressant response in MDD (Kishi et al, 2009c), despite its key role in circadian rhythm regulation. These findings deserve further investigation in larger samples considering homogeneous clinical phenotypes and circadian subphenotypes. Our results also corroborate previous reports suggesting that genetic variation in clock genes could lead to abnormal functionality of the circadian pacemaker underlying the biological and social rhythm disruptions, which have been implicated in triggering mood episodes in both unipolar and bipolar MD (Grandin et al, 2006).

In addition to the involvement of circadian genes in sleep homeostasis mentioned above, changes in the molecular clock could affect the circadian rhythms of secretion and function of other elements that have also been implicated in the pathophysiology of MD (Nestler et al, 2002). As circadian alterations have been described in monoaminergic secretion (Morin, 1999; Wirz-Justice, 1987), in hypothalamic–pituitary–adrenal axis hormones (Beck-Friis et al, 1985; Linkowski et al, 1987), and in immune system components such as cytokines (Humphreys et al, 2006), the chronobiological disturbances may act as an essential link in the pathogenesis of MD, complementing other etiopathogenical pathways and integrating its multidirectional interactions.

Certain methodological issues in this study deserve comment. First, population stratification may be a confounding factor in case–control genetic studies. We tested this possibility and found no evidence of stratification, suggesting that it does not introduce significant bias into our analyses. Second, genetic factors that contribute to MD probably comprise both common genetic variants with a small effect, such as the ones we have identified in this study, and rare sequence variants with a larger effect (Campbell and Manolio, 2007). Given that we did not include a mutation scan to assess rare variants in our study, we cannot rule out the participation of rare variants with functional effects in LD with the SNPs we identified for association in circadian genes with MD. It would therefore be necessary to resequence the coding and putative transcriptional regulatory genomic regions of CRY1, NPAS2, CLOCK, VIP, and VIPR2 in order to disentangle the genuine genetic effect sizes that contribute to MD phenotypes. Third, MD cases from clinical samples of tertiary sources may differ from community-based cases due to selection bias, and thus make comparisons with other studies difficult. Finally, we studied a set of 19 circadian clock genes, including the core clock genes and a large number of genes involved in their regulation and the molecular clock machinery (see Supplementary Figure S1 and Supplementary Table S1), but we cannot discard that other circadian or clock-controlled genes not studied in this study, such as SIRT1 (Asher et al, 2008) and PROKR2 genes (Prosser et al, 2007), could also be involved in the pathophisiology of MD. In accordance, a recent study conducted in Japanese population reports an association of PROKR2 gene with MD (Kishi et al, 2009e).

In conclusion, our data support the role of the circadian system in the genetic susceptibility to MD, pointing toward a contribution of circadian gene-specific effects to major MD polarity; some genetic variants in certain circadian genes seem to be implicated in depressive pole liability and others in manic pole liability. Our findings underline the importance of considering circadian system as a therapeutic target for the treatment of MD and developing specific therapeutic agents for unipolar and bipolar mood episodes. Moreover, the identification of specific genetic markers of vulnerability to bipolarity could have importance in the clinical management of major depression in the future, with prognostic significance and therapeutic implications. These results argue in favor of the recording of phenotype data related to circadian rhythms in MD that will allow the performance of additional analyses considering more refined clinical subphenotypes in our sample under the specific circadian hypothesis. Further studies in independent samples are warranted to confirm our present findings.