1. Cell Biology
Download icon

The genomic landscape of human cellular circadian variation points to a novel role for the signalosome

  1. Ludmila Gaspar
  2. Cedric Howald
  3. Konstantin Popadin
  4. Bert Maier
  5. Daniel Mauvoisin
  6. Ermanno Moriggi
  7. Maria Gutierrez-Arcelus
  8. Emilie Falconnet
  9. Christelle Borel
  10. Dieter Kunz
  11. Achim Kramer
  12. Frederic Gachon
  13. Emmanouil T Dermitzakis
  14. Stylianos E Antonarakis
  15. Steven A Brown  Is a corresponding author
  1. University of Zurich, Switzerland
  2. University of Geneva, Switzerland
  3. Laboratory of Chronobiology, Germany
  4. University of Lausanne, Switzerland
  5. Institute of Physiology, Charité-Universitätsmedizin Berlin, Working Group Sleep Research & Clinical Chronobiology, Germany
Research Article
  • Cited 3
  • Views 1,174
  • Annotations
Cite this article as: eLife 2017;6:e24994 doi: 10.7554/eLife.24994

Abstract

The importance of natural gene expression variation for human behavior is undisputed, but its impact on circadian physiology remains mostly unexplored. Using umbilical cord fibroblasts, we have determined by genome-wide association how common genetic variation impacts upon cellular circadian function. Gene set enrichment points to differences in protein catabolism as one major source of clock variation in humans. The two most significant alleles regulated expression of COPS7B, a subunit of the COP9 signalosome. We further show that the signalosome complex is imported into the nucleus in timed fashion to stabilize the essential circadian protein BMAL1, a novel mechanism to oppose its proteasome-mediated degradation. Thus, circadian clock properties depend in part upon a genetically-encoded competition between stabilizing and destabilizing forces, and genetic alterations in these mechanisms provide one explanation for human chronotype.

https://doi.org/10.7554/eLife.24994.001

Introduction

A biological ‘circadian’ clock governs nearly all aspects of human behavior and physiology in synchrony with the geophysical day. Although the basic mechanism of this clock is highly conserved across evolution, humans show widely varying phases of behavior (chronotypes) ranging from very early (so-called ‘larks’) to very late (‘owls’). In rare instances of familial circadian disorders such as advanced sleep phase syndrome, it has been possible to identify specific mutations in dedicated ‘clock proteins’ that affect human circadian behavior and the period of the underlying clock (Toh et al., 2001; Xu et al., 2005). More recently, these investigations have been complemented by genome-wide association (GWAS) studies using questionnaire-based methods to interrogate genotyped cohorts about diurnal preference. These strategies have unearthed significant polymorphisms in regions containing known clock genes as well as other loci, but linking novel loci to biological mechanisms remains challenging (Allebrandt and Roenneberg, 2008; Hu et al., 2016; Jones et al., 2016; Lane et al., 2016).

In cell- and tissue-based studies, it has been possible to link GWAS-based insights to biochemical mechanisms using eQTLs (expression quantitative trait loci) based upon transcript levels or other cellular parameters. Typically, because of the labor involved in collecting and analyzing cellular material, the cohorts used in such investigations are much smaller (hundreds of samples rather than >10,000 in conventional GWAS). Partially compensating for this deficiency is the relative simplicity of the analytical system and the precision of expression trait measurement (Nica and Dermitzakis, 2013). Such approaches have been used successfully in a variety of contexts, identifying susceptibility loci for autoimmune disease (Kochi, 2016), for HIV infection (Loeuillet et al., 2008), and multiple other diseases. Since individual common alleles typically have rather small effect sizes, eQTL analysis provides a cellular methodology to analyse disease-relevant molecular phenotypes in simple systems that can amplify their effects (Hou and Zhao, 2013).

Since circadian clock function is cell-autonomous, it would be ideally suited to eQTL-based investigations. The fundamental unit of circadian timekeeping is comprised of feedback loops of transcription and translation of dedicated ‘clock genes’, coupled to parallel and perhaps even independent cycles of post-translational modification. These highly conserved clocks are hierarchically organized in mammals including man, with a ‘master’ clock in the suprachiasmatic nuclei (SCN) of the hypothalamus driving diurnal behavior and physiology both directly and indirectly via ‘slave’ oscillators of similar mechanism elsewhere in the brain and body (Brown and Azzi, 2013).

In the recent past, circadian clocks in cultured fibroblast cells have been used as a workhorse for understanding circadian function. Importantly, interfering with these clocks produces clock phenotypes similar to those present in a corresponding whole-animal context (Brown et al., 2005). Questionnaires about human daily behavior are even predictive of fibroblast clock properties from the same individuals: in general, cells from ‘larks’ show a short period, and those from ‘owls’ a longer period (Brown et al., 2008). Herein, we have exploited these cellular proxies to explore the genetic basis of variation in human circadian clock function.

Results

eQTL GWAS identifies polymorphisms influencing cellular circadian function

To assay circadian clock function in cultured cells, differently-phased clocks in the cells of the culture dish are synchronized by rapid chemical induction of a clock gene, and then circadian gene expression is monitored for several days, typically via a synthetic reporter (Gaspar and Brown, 2015). With this aim, we transduced the 159 umbilical cord fibroblasts of the Genecord II library (Borel et al., 2011; Dimas et al., 2009) – each genotyped for 2.5 million common single-nucleotide polymorphisms, or SNPs -- with a lentiviral BMAL1-luciferase circadian reporter, then synchronized clocks with dexamethasone, and measured acute induction of PER1 and PER2 gene expression via qPCR, as well as the period, phase, and amplitude of subsequent circadian oscillation via real-time in vitro bioluminescence recording. For each parameter, large inter-individual differences were observed (Figure 1—figure supplement 1a,b), as has been reported previously (Brown et al., 2005). Normalized data were used as quantitative traits for genome-wide association, identifying for each clock parameter polymorphisms at different confidence levels up to p=10−7 (Figure 1a, see individual Manhattan and Q-Q plots in Figure 1—figure supplement 1c–g, Figure 1—figure supplement 2; all significant SNPs are listed in Figure 1—source data 1).

Figure 1 with 3 supplements see all
eQTLs influencing circadian function.

(a) Manhattan plot of identified polymorphisms for all measured traits. Blue line, threshold for suggestive candidates, p<10−5, approx. FDR ≤ 0.1. Genes likely associated with the most significant alleles are shown in red. (b) For all SNPs at the indicated p-value (x-axis), Tukey Boxplot of the distribution of difference in allele frequency between self-declared ‘larks’ and ‘owls’ (y-axis, absolute value of (frequency in larks – frequency in owls); x-axis, p-value threshold used for the test; number of alleles at each stringency indicated underneath). (c) Equivalent chart for Tag-SNPs only, that is one SNP per hapblock. See also Figure 1—figure supplement 1 for distributions and Manhattan plots for individual traits; see Figure 1—figure supplement 2 for quantile-quantile plots associated with each Manhattan plot; see Figure 1—figure supplement 3 for density plots of allele frequency distributions in extreme chronotypes. Figure 1—source data 1 lists all alleles identified with p<10−5 for each trait, as well as the closest gene.

https://doi.org/10.7554/eLife.24994.002
Figure 1—source data 1

SNPs associated with (a) PER1 and (b) PER2 expression, (c) Amplitude, (d) Phase, and (e) Period.

The closest annotated gene is indicated for each SNP, as well as the unadjusted GWAS p-value.

https://doi.org/10.7554/eLife.24994.006

We next compared SNPs across the traits that we identified. Considering alleles achieving ‘suggestive’ genome-wide significance (p<10−5), we see correlation between multiple pairs of traits, notably PER1 and PER2 expression (p=0.024), period and phase (p=1.48e-10), and period and amplitude (p=0.009) (Supplementary file 1). Given the close homology and overlapping function of the PER proteins, and the known correlation between period length and phase, these correlations are expected. Surprisingly, however, the most significant alleles in any one category are not among the most significant alleles in another, suggesting relatively independent genetic regulation of key circadian clock parameters (‘state variables’) for small changes, and the resilience of the overall circadian mechanism to small perturbations.

While multiple suggestive alleles were identified, it should be noted that none achieved ‘genome wide’ significance (often reckoned at p<5×10−8), as might be expected from a cell-based study of moderate sample size and 2.5 million SNPs tested. Therefore, we undertook further experiments to demonstrate the relevance of our results.

Putative significant alleles are enriched in extreme chronotypes

To explore the significance of our results for human behavior and to probe the confidence levels required to elicit meaningful data, we examined the enrichment of SNP alleles in individuals of extreme chronotype after genotyping a previously-characterized cohort of 35 ‘larks’ and ‘owls” (Brown et al., 2008). We reasoned that for SNPs identified at increasing stringency of p-value, those affecting behavior would show an allelic frequency that varied between larks and owls. Consistent with this hypothesis, allele enrichment analysis shows that differential allelic frequency among SNPs identified at relaxed confidence levels up to p>10−4 is very low, and then increases at p<10−5 or p<10−6, plotted either for all significant SNPs (Figure 1b), or as single tag-SNPs per hapblock (Figure 1c). The same data is represented as a density distribution in Figure 1—figure supplement 3. Here, it can be seen that whereas most alleles have a differential frequency of 0 between larks and owls, this rises progressively to 0.1 for the 22 alleles identified as significant at a threshold of p<10−6.

Genes associated with positive SNPs affect circadian clock function

As a next step, for each SNP with FDR < 0.1, (using multiple-comparison correction for the 2.5 million SNPs tested) we identified putative genes whose transcription start site lay within 100 kb of the indicated SNP (76 genes). We chose this threshold to correspond to ‘suggestive’ significance (uncorrected p<1×10−5) in genome-wide association: since it is estimated that only a tenth of the SNPs on the high-density chip that we used would be necessary to capture all common variation in individuals of European origin (International HapMap Consortium, 2005), a corrected FDR is therefore up to tenfold ‘over-corrected’. We therefore preferred to select genes at a relatively relaxed FDR q-value, and subsequently apply further criteria.

For genes more than 50 kB away from a significant SNP, the further criterion of significant association of the SNP region with the gene of interest via chromatin conformation capture (3C) was applied, using the HiView fibroblast dataset (Xu et al., 2016). In total, a list of 59 genes putatively affecting human chronotype was obtained in this way. For example, the two most significant SNPs uncovered in our screen, rs920400 (p=1.0E−7) and rs10195385 (p=1.9E−7) affecting the expression of PER2, lay in the conserved enhancer region of the COPS7B gene, encoding a subunit of the human COP9 signalosome (Fang et al., 2008). These SNPs show a strong 3C signal with the COPS7B gene region (Figure 2a), supporting their annotation as affecting COPS7B expression. (Please see later results for more detailed biochemical and genetic analysis of the effects of these two alleles.)

Figure 2 with 1 supplement see all
Assessment of genes associated with identified eQTLs.

(a) Genome-wide high-throughput chromosome conformation capture (Hi-C) for probes spanning rs920400 (blue) and rs10195385 (black). Y-axis, read counts. X-axis, genomic location of reads. Gene neighborhood is shown in the inset. (b) RSA-score distribution of period lengths of U2OS cell cultures transfected with RNAi hairpins targeting each gene in the genome (grey) or in significant genelist (blue). Y-axis, relative density (area sums to 1); X-axis, log RSA score. Left of 0, period length less than normal; right of 0, period length greater than normal. Blue circles, co-plotted scatterplot of individual values for significant genelist; black circles, corresponding sample set drawn randomly from the whole genome. (c) RSA-score density distribution of amplitude. (X and Y axes as in (b).) An annotated Significant Genelist is given in Supplementary file 1. See also Figure 2—figure supplement 1 for circadian clock properties of cells in which significant genes are targeted by RNAi, as well as representative raw data and cumulative density plots.

https://doi.org/10.7554/eLife.24994.007

A second way of looking for associations between SNPs and genes relies upon correlating SNPs with changes in transcription (Regulatory Trait Concordance). Since transcriptomic data from these cells were already available (Gutierrez-Arcelus et al., 2015), we used it in order to also search for potential cis- and trans-effects via Regulatory Trait Concordance (Nica et al., 2010). This procedure uncovered one additional gene (PPM1B) at these significance criteria. (All genes and their corresponding SNPs are listed in Supplementary file 2).

Surprisingly, none of the significant SNPs were present in ‘canonical’ clock genes of the circadian transcription-translation feedback loop. Therefore, to verify the significance of our findings, we conducted RNA interference (RNAi) targeting each of the genes that we identified using the human model circadian cell line U2OS, and looked for changes in circadian period and amplitude in vitro. Of these 60 genes, depletion of 28 genes showed a significant difference in circadian period of at least one standard deviation from the Z-scored mean, and 29 genes showed similar alterations in circadian amplitude. (Data for all genes tested and representative circadian profiles are shown in Figure 2—figure supplement 1a,b). To quantify the significance of these results, we used Redundant siRNA Activity (RSA) scores to attribute an overall p-value to siRNAs tested for each gene (König et al., 2007). Thus, for each clock trait – period and amplitude – we obtained two distributions: p-value probability of increase, and p-value probability of decrease of these traits for each gene identified in our screen. We then compared these to probabilities obtained for the transcriptome as a whole. Overall, the distribution of p-values for period obtained from hairpins targeting genes identified in our screen shifts outward (toward greater significance) and broadens, that is these genes have an increased probability of affecting the circadian clock compared to random genes (Figure 2b, Kruskal-Wallis Test p=0.14 reduced and p=0.0016 increased period; the combined cumulative distribution is shown in Figure 2—figure supplement 1c, Kolmogorov-Smirnov Test p=0.07). Analogous effects, but stronger, are observed for amplitude (Figure 2c, Figure 2—figure supplement 1c, p=0.0002 reduced and p=0.034 increased, p=0.0013 combined). Interestingly, a skew to the left is also visible though not significant (Figure 2—figure supplement 1c), perhaps reflecting the greater probability of reducing circadian amplitude by depleting identified genes than increasing it.

Positive SNPs highlight protein catabolism as essential to human circadian variation

To try to understand the broad biochemical pathways underlying the genes that were identified, we first applied gene ontology enrichment analysis (Wang et al., 2013) to the list of 60 genes that we obtained in Supplementary file 2. The results, depicted as an abbreviated Gene Ontology diagram in Figure 3a (full diagram shown in Figure 3—figure supplement 1; individual genes and categories listed in Figure 3—source data 1), highlight in particular signaling, development, and macromolecular metabolism as significant, with half of all positive genes involved in metabolic processes (Figure 3b). To independently confirm these relationships and explore them in more detail, we next used the completely independent approach of a gene set enrichment algorithm (GSEA) adapted for GWAS studies, iGSEA4GWAS (Zhang et al., 2010). This algorithm begins by attributing genes to positive SNPs using a distance-based criterion (the same 100 kb that we used previously), looks for pathways with high proportions of significant genes at relaxed stringency criteria (p<10−4), and imputes a global significance to this association weighted by SNP p-values of member genes. Applying this method to our SNP data, several pathways of potential significance emerge, including Hedgehog signaling, Cellular Oxidoreductase Activity, Cell Cycle G1/S, and Transcriptional Repressor Activity. A list of these genes and SNPs is presented as Supplementary file 3. (A list of categories, all genes in each, and significant genes in each is available in Figure 3—source data 2). Five of the top six identified pathways are involved in macromolecular catabolism. The full list of pathway annotations from GSEA analysis is shown in Figure 3—figure supplement 2a. We also show an abbreviated list, eliminating redundancies in pathway annotations, in Figure 3c.

Figure 3 with 2 supplements see all
Positive SNPs highlight protein catabolism as essential for human circadian variation.

(a) Abbreviated Gene Ontology diagram from the Gene SeT Analysis Toolkit (WebGestalt), showing biological process categories containing genes (Supplementary file 1) associated with positive SNPs. Boxes with red statistical description represent enriched categories, with the number of associated genes and corresponding multiple comparison-adjusted p-value. (Individual pathways and associated genes are annotated in Figure 3—source data 1.) (b) Biological processes categorization of all identified genes (GOSlim tool). Y axis, number of genes in the category. (c) The most significant nonredundant pathways associated with positive alleles, identified by GeneSet Enrichment Analysis (GSEA4GWAS toolkit; a full chart including closely-related pathways is given in Figure 3—figure supplement 1; All genes are listed in Supplementary file 2; individual pathways and a clickable list of associated genes and SNPs is given in Figure 3—source data 2.) (d) Manhattan plot of all alleles p>10−3. GSEA-uncovered genes related to protein catabolism are indicated, and associated polymorphisms are circled. (e) RSA-score distribution of period lengths of U2OS cell cultures transfected with RNAi hairpins targeting each gene in the genome (grey), or targeting genes in (d) (blue). Circles, co-plotted scatterplot of individual values. X and Y axes as in Figure 2b,c. (f) RSA score distribution for circadian amplitude observed in U2OS cell cultures transfected with RNAi hairpins targeting each gene in the genome (grey), or targeting genes in (d) (blue). Circles, co-plotted scatterplot of individual values. X and Y axes as in Figure 2b,c. See also Figure 3—figure supplement 1 for a full list of GSEA-identified pathways, as well as a list of RNAi-targeted genes associated with protein catabolism, and the circadian clock properties of the cells in which this targeting occurred.

https://doi.org/10.7554/eLife.24994.009
Figure 3—source data 1

Multiple-comparison-adjusted p-values and significant associated genes for each pathway identified in Figure 3a (and in more detail in Figure 3—figure supplement 1).

https://doi.org/10.7554/eLife.24994.012
Figure 3—source data 2

Multiple-comparison-adjusted p-values, false discovery rates, and significant associated genes for each pathway identified in Figure 3—figure supplement 2a.

Significantly associated genes, as well as all genes in the pathway not significantly associated, can be perused via the ‘view detail’ links for each pathway.

https://doi.org/10.7554/eLife.24994.013

Since by far the most significant identified category (5 out of the top 6 terms) was protein catabolism, we chose to explore further the significance of alleles in this group. All alleles and their affected genes from this category are shown on a Manhattan plot in Figure 3d, and listed in Figure 3—figure supplement 2b. The RSA score distribution of these genes showed an increased probability of long- and short-period phenotypes – seen as an outward-shifted and broadened distribution (p=0.009 overall, 2-way Kolmogorov-Smirnov Test; p=0.044 shorter, p=0.0008 longer, Kruskal-Wallis Test; Figure 3e, Figure 3—figure supplement 2d), and a significant skew toward reduced circadian amplitude (p=0.024 reduced, p=0.47 increased, Kruskal-Wallis Test, Figure 3f, more easily seen as a cumulative distribution in Figure 3—figure supplement 2e). Thus, the results of our genetic screen support the hypothesis that a major reason for human circadian variation at a cellular level relates to global alterations in protein stability.

COPS7B influences PER2 expression and clock function.

Consistent with this idea, the top two SNPs identified in our screen are in putative regulatory regions of the COPS7B gene (as well as 24 other alleles in the same hapblock with p<10−6). This protein is part of the COP9 signalosome, which has been implicated both in regulating ubiquitin-mediated proteolysis (Pick and Bramasole, 2014), and in regulating circadian behavior of Neurospora crassa and Drosophila melanogaster (He et al., 2005; Knowles et al., 2009). Therefore, understanding how COPS7B affects circadian function could provide important clues to the mechanisms of genetic influence upon human circadian clocks.

To access the effect of COPS7B allelic identity upon expression of a clock gene, we independently measured PER2 expression values as a function of rs920400 genotype, and found that average PER2 expression in the rs920400 AA genotypes were significantly higher (p=0.00026) than the GG genotype (Figure 4a).

Figure 4 with 1 supplement see all
COPS7B influences PER2 expression and cellular clock function.

(a) Expression levels of PER2 stratified by rs920400 genotype. (AA vs. GG genotype, T-test p=0.00026). Insets - Correlation between PER2 and COPS7B gene expression within each genotype. (b) Reduction of PER2 expression (grey) and COPS7B expression (black) in HEK293T and U2OS cells transfected with three short interfering RNAs targeting COPS7B (c) Raw bioluminescence profiles of U2OS cells cotransfected with a Bmal-luciferase circadian reporter and siRNAs targeting COPS7B (dotted lines) or scrambled siRNA (solid lines). (d) Reduction in circadian amplitude by short interfering RNAs targeting COPS7B (p=0.0021, student T-test). These siRNAs also lengthen the circadian period by 1.15 hr (p=3.27094E−05). See also Figure 4—figure supplement 1 for similar RNAi assays performed upon the COPS4 gene.

https://doi.org/10.7554/eLife.24994.014

As mentioned previously, this SNP falls in a conserved putative enhancer region for its closest gene, COPS7B, lying 30 kb away. To demonstrate a potentially functional relationship between COPS7B and PER2, for each allele we plotted the relationship between COPS7B expression and PER2 expression, and demonstrated strong positive correlations within each genotype (AA: r = 0.9327, p=0.0007; AG: r = 0.9455, p=0.0004, GG: r = 0.8666, p=0.0054) (Figure 4a, insets).

It should be noted that here and in the screen itself, a single timepoint was used to investigate PER2 expression, corresponding to its expected first peak after cell line synchronization. Thus for any single subject, absolute changes in level could also reflect changes in timing. Therefore, we next explored the origins of this effect more fully using exogenous assays. We transfected the circadian model cell line U2OS:Bmal1-luc (Maier et al., 2009) with three different short interfering RNAs targeting COPS7B, reducing its expression by 68%. These siRNAs correspondingly reduced PER2 expression (Figure 4b), and also decreased circadian amplitude more than 10-fold, p=2.93×10−7 and lengthened circadian period by 1.15 hr, p=3.27×10−5 (Figure 4c,d), even as cells continued to divide viably. Showing that the modulation of circadian clock function is achieved by the COP9 signalosome rather than by COPS7B independently, a an even more severe phenotype was achieved by targeting COPS4, an adjacent COP9 signalosome subunit: circadian amplitude was reduced more than 10-fold in the first day, p=3.59E−7 and period initially lengthened by 8.8 hr (p=4.02E−5), followed by arrhythmicity, (Figure 4—figure supplement 1a,b). Because this phenotype is the opposite of the shorter period observed by eliminating either PER protein directly (Zheng et al., 2001), it also led us to suspect that the effects upon PER2 transcription were indirect.

COPS7B interacts physically with clock proteins and stabilizes BMAL1

To demonstrate physical interactions of the COP9 signalosome with the circadian clock machinery, we transfected HEK 293 T cells – which do not possess a functional circadian clock, thereby eliminating the variable of time – with constructs expressing FLAG epitope-tagged clock proteins BMAL1, PER2, and CRY1. Subsequently, each lysate was immunoprecipitated with either anti-FLAG or anti-COPS7B antibody, and then Western-blotted and probed with the opposite. Interactions between COPS7B and all three clock proteins were observed bidirectionally (Figure 5a,b), while no corresponding clock proteins were immunoprecipitated by FLAG antibody in the absence of transfection (Figure 5—figure supplement 1).

Figure 5 with 2 supplements see all
COPS7B stabilizes the essential clock protein BMAL1.

(a–b) Coimmunoprecipitation assay of HEK293T cells transfected with plasmids expressing FLAG epitope-tagged clock proteins (BMAL1-flag, PER2-flag and CRY1-flag). *Unspecific background bands. (a) immunoprecipitated with antibodies against FLAG or (b) against endogenous COPS7B. (c) Pulse-chase analysis of BMAL1-flag protein stability in absence (left; black diamonds in graph below) or presence (right; grey squares in graph below) of siRNAs targeting COPS7B in HEK293T cells. All cells were incubated with 35S-labelled methionine-cysteine for 1 hr, and chased with excess unlabeled cysteine for 0 hr, 2 hr, 4 hr, or 6 hr before immunoprecipitation. Upper panel, representative radioblot; lower panel, quantification from 3 experiments ± s.d., expressed as percentage of labeled immunoprecipitated protein (relative to 0 hr) at indicated time. *p=0.04999, Student T-test. (d) Immunoprecipitation experiments between BMAL1 and COPS7B performed in U2OS cells. See also Figure 5—figure supplement 1 for negative control experiments upon untransfected cells; see Figure 5—figure supplement 2 for pulse-chase analysis of PER2-flag and CRY1-flag protein stability in the presence or absence of siRNAs targeting COPS7B.

https://doi.org/10.7554/eLife.24994.016

The COP9 signalosome has been implicated previously in antagonizing ubiquitin-mediated degradation of proteins (Korczeniewska and Barnes, 2013; Fernandez-Sanchez et al., 2010). Therefore, we next measured the half-life of transfected FLAG-tagged circadian clock proteins by metabolic pulse-chase protein labeling in the presence and absence of siRNAs targeting COPS7B. The timecourses of degradation of CRY1 and PER2 remained unchanged (Figure 5—figure supplement 2), but the half-life of BMAL1 was reduced to 4.4 hr when compared to control whose half-life was 7.3 hr. (Figure 5c). Therefore, we hypothesized that the bona-fide interaction between COP9 and the circadian clockwork was with BMAL1, even if biochemical interactions with other clock proteins could be found in transfection-based assays..

COPS7B and the COP9 signalosome are rhythmically imported into the nucleus

Having demonstrated an effect of COP9 signalosome impairment upon BMAL1 protein stability in HEK 293 T cells, we next examined the same phenomenon in U2OS cells, which possess robust circadian rhythmicity. In these cells, interaction was observed with BMAL1 only (Figure 5d). A subsequent circadian timecourse using synchronized cells harvested at different times of day illustrated one possible cause for the selectivity of this interaction: by Western blot, COPS7B shows circadian nuclear abundance in phase with BMAL1 but antiphase to PER and CRY proteins (Figure 6a). It should be noted that this result does not completely explain why interactions with other clock proteins would not be observed when they are exogenously transfected, and suggests therefore that time-specific cofactors might be required for the interaction.

Figure 6 with 1 supplement see all
The signalosome is rhythmically imported into the nucleus.

(a) BMAL1, PER2 and COPS7B RNA levels in U2OS cells at different times of day, as determined via qPCR. Inset - Western blot analysis of nuclear COPS7B protein under identical conditions. Top panel, representative raw blot; bottom panel, quantification. (b) RNAseq transcript levels of genes encoding COP9 signalosome subunits in mouse liver at different times of day. (c) Protein levels of COP9 signalosome subunits measured by MS in mouse liver nuclei at different times of day. Red line, protein levels of circadian clock protein BMAL1. See also Figure 6—figure supplement 1 for ms-based relative abundance analysis of different COP9 signalosome subunits, as well as protein levels of proteasome subunits measured by MS in mouse liver nuclei at different times of day.

https://doi.org/10.7554/eLife.24994.019

To demonstrate the relevance of our findings in vivo, we turned to a mouse liver model. Examining a mouse liver transcriptomic dataset (Atger et al., 2015), the RNAs encoding most COP9 signalosome subunits showed arrhythmic transcription, but Cops7b transcripts were markedly circadian (Figure 6b), and also expressed at the lowest level (Figure 6—figure supplement 1a). Since our COPS7B antibodies were specific to human, we used isotope-labeled mass spectrometry to examine nuclear accumulation of COP9 signalosome components. All COP9 subunits are imported rhythmically into liver nuclei in phase with the peak of BMAL1 protein accumulation, and slightly in advance of it (Figure 6c), implying a circadian role for COP9 at this time of day. Thus, both rodent and cellular studies support the idea that changes in COP9 expression could alter circadian function.

Discussion

A salient feature of human circadian behavior is its difference across individuals.

Previous studies have focused upon the importance of ‘canonical’ clock genes – core members of the circadian transcription-translation feedback loop – in these variations. Specific clock gene mutations been identified both in families (Toh et al., 2001; Xu et al., 2005) and in the broader population (Allebrandt and Roenneberg, 2008; von Schantz, 2008) that could affect human behavior. Moreover, three recently-published large questionnaire-based GWAS studies of circadian behavior have also highlighted polymorphisms in the region of known clock genes as a sizeable portion of the signals they observe (Hu et al., 2016; Lane et al., 2016; Jones, 2015). Taking a cellular rather than a behavior-based GWAS approach to circadian variation, we find evidence of a wider biochemical truth: a major source for differences in clock function at a cellular level is in fact based upon the broad mechanism of protein catabolism. Surprisingly, there is no overlap between top candidate alleles in their studies and ours. Moreover, among these three studies, so far only meta-analysis of grouped data has been performed, and no evidence was seen for other previously implicated alleles. Thus, the possible space for novel discovery remains large. Moreover, it is unclear to what extent any given QTL affecting cellular clock properties such as those that we investigate would propagate directly to effects upon human behavior.

As a feedback loop, the circadian clock is critically dependent upon protein stability to time its oscillations. Therefore, our conclusions make mechanistic sense. Indeed, both genetic linkage studies cited above (Toh et al., 2001; Xu et al., 2005) have highlighted the importance of post-translational modification of clock proteins to regulate ubiquitin-proteasome pathway-dependent protein stability. RNAi screens (Maier et al., 2009; Sathyanarayanan et al., 2008; Zhang et al., 2009) and mouse mutant analyses (Godinho et al., 2007; Siepka et al., 2007) have also emphasized the global importance of proteasome-mediated degradation to circadian function. Here, by using human genetics as a tool for mechanistic discovery, we have demonstrated a novel role for the COP9 signalosome in the stability of the essential circadian clock protein BMAL1, a transcriptional activator of the PER2 gene (Gekakis et al., 1998) and a critical player in the circadian process. In fact, the entire signalosome is imported into the nucleus both in vitro and in vivo in time with the maximal abundance of BMAL1, enhancing its stability. Thus, the primary effect is a competition between signalosome and proteasome to stabilize or destabilize BMAL1, the former in time-specific fashion. Although we could detect interactions with multiple transfected clock proteins in nonrhythmic HEK cells, in rhythmic U2OS cells only BMAL1 interaction was detectable, suggesting further that other time-of-day-specific factors might be required for this interaction.

For the ‘negative limb’ of the circadian clock, such a phenomenon has already proven essential to the regulation of circadian period. Competition between two different E3 ubiquitin ligases, FBXL3 and FBXL21, respectively destabilize and stabilize CRY proteins (Hirano et al., 2013; Yoo et al., 2013). Analogous effects can be observed for PER2 protein, this time via competition between different phosphorylation events (Vanselow et al., 2006; Xu et al., 2007). Competition for phosphorylation sites and ubiquitin ligases is also a key determinant of the speed of the Drosophila circadian oscillator (Chiu et al., 2011; Chiu et al., 2008). With these genetic studies, we have unearthed evidence of a conceptually similar but mechanistically different regulation of the ‘positive limb’ via BMAL1 stability.

It is worth noting that although knockdown of COP9 subunits altered circadian period, phase, and PER2 expression, the COPS7B-associated alleles rs920400 and rs10195385 that we genetically identified were associated significantly only with PER2 expression in vitro: rs920400/PER2: p=1.0677E−7, rs10195385/PER2: p=1.9463E−7. For the same SNPs, association with other clock properties was not significant -- rs920400/PER1: p=0.202572306, rs10195385/PER1: p=0.198630802; rs920400/AMPLITUDE: p=0.791990033, rs10195385/AMPLITUDE: p=0.944165947; rs920400/PHASE: p=0.489070695, rs10195385/PHASE: p=0.522821404; rs920400/PERIOD: p=0.602715502, rs10195385/PERIOD: p=0.700789177. More generally, although there was detectable and significant genotypic correlation between period and phase, and PER1 and PER2 expression in our screen (Supplementary file 1Figure 1—figure supplement 3), genetic correlation among these different phenotypes was nevertheless weaker than we expected. We attribute this weakness to the idea that small changes in cellular clock properties can be robustly tolerated without affecting other clock parameters. Such an idea has also been supported by cellular dosage studies individual clock proteins (Baggs et al., 2009).

The COP9 signalosome has been shown previously to inactivate many cullin-RING-ubiquitin E3 ligases, probably by removing ubiquitin-like NEDD8 from cullins (so-called ‘de-neddylation’) (Dubiel et al., 2015; Cope et al., 2002). It has also been demonstrated to interact with casein kinase 2 (Uhle et al., 2003), which itself has a circadian role in mammals (Maier et al., 2009). Since the proteasome itself is constitutively present in the nucleus (Figure 6—figure supplement 1b) but the signalosome is circadian, it could thereby provide a more general mechanism for circadian regulation of nuclear protein stability. For example, in other eukaryotes, it has been suggested that the COP9 signalosome plays a role in the ubiquitin-mediated degradation of the clock protein TIMELESS in Drosophila melanogaster (Knowles et al., 2009), and Frequency in Neurospora crassa (He et al., 2005). Notably, the role of COP9 in Drosophila is also stabilizing, but opposite in effect: by stabilizing the ubiquitin ligase itself, COP9 thereby promotes degradation of TIMELESS.

Mechanisms related to protein stability are the most prominent identified here -- not only the signalosome, but also kinases, phosphatases, ubiquitin ligases, and regulatory proteins. However, it was not the only one. For example, the Hedgehog signaling pathway is represented at high significance, not only through casein kinase 1 homologs – a protein family known to have circadian function (Cheong and Virshup, 2011) – but also through WNT secreted signaling molecules. These were recently shown to play a role in cell cycle-circadian clock synchronization in intestinal organoids (Matsu-Ura et al., 2016). Cell cycle proteins are also represented – logically, since the cell division and circadian clock cycles are robustly coupled (Bieler et al., 2014), but nevertheless mechanistically not well understood. Likewise, cellular redox has come to recent attention in the context of circadian clocks (Putker and O'Neill, 2016), but specific mechanisms remain to be defined. Overall, human genetic studies have considerable power to discover novel physiologically relevant circadian mechanisms, and cell-based approaches make this discovery process even easier. In the case of circadian biology, such investigations could both provide new targets directly relevant to human health, and help to explain clock connections to disease.

Materials and methods

Ethical approval

Request a detailed protocol

All human samples used in these studies were obtained after approval of all protocols and procedures by the relevant responsible authorities (of the University Hospital Geneva, CH; and Charité Universitätsmedezin, Berlin, DE), and prior written informed consent was obtained from all subjects.

Tissue isolation, cell culture and genotyping

Request a detailed protocol

Umbilical cord fibroblasts were collected from 159 newborns of Western European origin as described previously (Dimas et al., 2009). Briefly, under sterile conditions, the cord tissue was finely cut in 1 ml DMEM containing 10% FCS, 1% antibiotics (Amimed, Basel, Switzerland), transferred to a T25 flask and cultured upside-down for 12 hr to allow cells to attach to the surface of the flask. Flasks were then turned and kept for about 1 week until fibroblast clusters appeared. Subsequently, fibroblasts were expanded with standard procedures. Similarly, fibroblasts from Western European adults of extreme chronotype (17 ‘owls’ and 11 ‘larks’) were cultured as described previously (Brown et al., 2008).

All fibroblasts were genotyped with the Illumina Omni2.5–8 and the Omni2.5Exome-8 chips. Population stratification was done by principal component analysis with phase 1 1000 genomes variants. The phasing was performed with the SHAPEIT v.2.r790 software (Delaneau et al., 2011) and the imputation with the IMPUTE2 v2.3.1 software using the phase 1 1000 genomes genotypes as a reference panel. Variants were filtered out according to these criteria: MAF <5%, Hardy-Weinberg probability <1e-6, call rate <97.5% and imputation score <0.4. A total of 5210911 variants were left after these filtering steps.

Fibroblast culture, viral transduction, synchronization, and measurement of circadian rhythms parameters

Request a detailed protocol

Cultivation of fibroblast cells, virus preparation, transduction, and selection of stable transformants were performed as previously described (Brown et al., 2005). Before the measurement, circadian rhythms in identically grown plates were synchronized with the synthetic glucocorticoid dexamethasone (Balsalobre et al., 2000). After washing twice with 1x PBS, cells were incubated in DMEM medium without phenol red, supplemented with 10% FBS, 0.1% gentamycin and 0.1 mM luciferin. Light emission was measured in a lumicycle photomultiplier device at 1 min intervals at 37°C, 5%CO2 for 5 days.

RNA preparation and gene expression quantification

Request a detailed protocol

Total RNA was extracted by using High Pure RNA Isolation Kit (Roche) following the instructions of the manufacturer. 500 ng of isolated RNA was transcribed with SuperScript II (Invitrogen, Carlsbad, CA) using random hexamer primers according to the manufacturer’s instruction. For quantitative PCR, 20 ng of cDNA was used and single transcript levels of genes were detected either by TaqMan probes used with the TaqMan PCR mix protocol (Roche, Basel, Switzerland) or with HOT FIREPol EvaGreen qPCR Mix Plus (Solice Biodyne, Tartu, Estonia) and an AB7900 thermocycler (Applied Biosystems, Foster City, CA). Primers used for detection of transcripts were as follows: BMAL1 forward: 5’-gaagacaacgaaccagacaatgag-3’, BMAL1 reverse: 5’- acatgagaatgcagtcgtccaa-3’, BMAL1 TaqMan probe:5’-FAM-tgtaacctcagctgcctcgtcgca-TAMRA-3’; PER1 forward: 5’-cgcctaaccccgtatgtga-3’, PER1 reverse: 5’-cgcgtagtgaaaatcctcttgtc-3’, PER1 TaqMan probe: 5’-FAM-cgcatccattcgggttacgaagctc-TAMRA-3’; PER2 forward: 5’-gggcagcctttcgactattct-3’, PER2 reverse: 5’ gctggtgtccaacgtgatgtact-3’, PER2 TaqMan probe:5’-FAM-cattcggtttcgcgcccggg-TAMRA-3’; GAPDH forward: 5’-cacatggcctccaaggagtaa-3’, GAPDH reverse:5’-tgagggtctctctcttcctcttgt-3’, GAPDH TaqMan probe:5’-FAM-tggaccaccagccccagcaaga-TAMRA-3’; BLASTICIDIN forward: 5’- gcgacggccgcatct-3’, BLASTICIDIN reverse: 5’-acaaggtcccccagtaaaatg-3’,

BLASTICIDIN TaqMan Probe: 5’-FAM-cactggtgtcaatgtat-TAMRA-3’ (FAM is 6-carboxyfluorescein, and TAMRA is 6-carboxytetramethylrhodamine);

COPS7B forward: 5’-ggctgggattagggtggttc-3’, COPS7B reverse:5’-atgtgggtacaggccttcctc-3’.

Data analysis of bioluminescence measurements and transcript levels

Request a detailed protocol

The circadian period length, phase, and amplitude were analyzed via the Waveclock R package (Price et al., 2008). To normalize circadian amplitude, the relative expression levels of BLASTICIDIN determined by qPCR for each subject were used. Data sets for the period length, phase, and amplitude were measured and plotted as average ±standard error from four replicates in three independent measurements for each subject. Peak expression levels of PER1 and PER2 likewise the levels of COPS7B and BLASTICIDIN transcripts were measured in four replicates for each tested subject.

Phenotype normalization

Request a detailed protocol

Prior to GWAS analysis, all data was corrected to eliminate systemic differences related to technical variability from cell passage and date of experiment. First, we averaged data coming from technical replicate of the same biological replicate. Second, we noted that some of our phenotypes (circadian period length, phase and amplitude) varied slightly but globally with the passage number of cell lines as well as with the date of the experiment. In order to eliminate effect of these confounding parameters we residualized our phenotypes using backward multiple regression model (lme function in nlme R package): starting from all variables included into the model we step by step eliminated less significant ones until all variables were nominally significant (p<0.05). Passage number we treated as ordinary numerical variables while date of experiments – as factors (dummy variables). Because all biological replications are non-independent, that is they belong to the same cell line, we used the cell line ID as a grouping variable in the lme function. Third, averaging residualized data coming from each biological replication we obtained the final phenotype value for each cell line. These values were used for all downstream analyses.

Genome-wide association study and statistical analysis

Request a detailed protocol

Independent GWAS were performed on the following normalized phenotypes: PER1 (140 samples), PER2 (141 samples), amplitude (142 samples), phase (147 samples) and period (147 samples). The GWAS were performed by using Spearman correlation between the genotypes and the phenotypes (continuous variable). All nominal p-values were reported on the Manhattan plots.

Regulatory Trait Concordance (RTC)

Request a detailed protocol

In order to link our GWAS hits with potential causal variants affecting the expression of genes, we used the RTC method described by Nica et al (Nica et al., 2010). GWAS variants with p-values smaller than 1E−4 were queried by the RTC method on a list of expressed quantitative trait loci (eQTL) found on GENCORD fibroblast samples (Gutierrez-Arcelus et al., 2015). The GWAS variants rs1878511 (p-value: 1E-4.089) and rs34041043 (p-value: 1E-4.089) had a RTC score bigger than 0.9 indicating a colocalization of the GWAS signals with the eQTL gene: PPM1B (SNP: rs4953137).

Allele enrichment analysis

Request a detailed protocol

This method was developed to show that alleles identified as affecting circadian clock properties in vitro would also affect human behavior in vivo. For each genotyped SNP, an allele frequency difference was calculated as the absolute value of (allele frequency in larks – allele frequency in owls). Larks and owls here indicate a group of individuals of extreme early and late chronotype identified in a previous study from our laboratory (Brown et al., 2008). The average frequency difference was then calculated for each group of alleles identified as significant in our genetic screen for a given p-value, and graphed either as a bar-and-whiskers plot (Figure 1b,c) or as a density plot (Figure 1—figure supplement 3).

RNA Interference assays

Request a detailed protocol

Lentiviral vectors expressing RNAi hairpins were purchased from Open Biosystems. Viral production, infection of U-2 OS cells, and subsequent measurement of circadian oscillations was performed as described in (Maier et al., 2009), n = 1/construct. Data analyses were also performed as described therein. Z-scored values for period and amplitude were calculated from individual RNAi constructs, and multiple RNAi constructs for a given gene were then concatenated using Redundant SiRNA Activity (RSA) scores calculated as described (König et al., 2007). Density distributions of RSA scores were calculated for a given genelist using all hairpins available in the Open Biosystems library, and compared to that for the entire transcriptome present therein using density, Kolmogorov-Smirnov, and Kruskal-Wallis tests from the R stats package. Density distributions were plotted either as mirror-image half-graphs showing probabilities of increase and decrease on the same axes, or as cumulative distributions of the combined set of both probabilities. (Practically, this was achieved by inverting the sign of the P value distributions reflecting increased period or amplitude.) For display purposes, the entire set of genelist data was shown as open circles, and compared to a random sample of 50 genes randomly chosen from the increased value distribution and 50 from the decreased value distribution. Note that RSA scores for hairpins showing different results for the same gene (e.g. one with increased and one with decreased period) will display a probability in both distributions, so the number of points observed can be greater than the total number of genes.

Bioinformatics

Chromatin conformation capture data was obtained from the HiView dataset and web-based search engine (Xu et al., 2016), and visualized on the UCSC genome browser http://genome.ucsc.edu/. Pathway analysis was done using the Webgestalt query tool (Wang et al., 2013) with default settings. GSEA4GWAS (Zhang et al., 2010) was used with 100 kB windowing and canonical pathways.

Cell lines

Request a detailed protocol

HEK293T and U2OS cells used for these studies were purchased from Sigma-Aldrich, St. Louis, MO. U2OS:Bmal1-luc cells were constructed as described (Maier et al., 2009). All lines were verified to be mycoplasma-free by PCR-based testing shortly prior to these experiments.

Transient transfection protocols

Request a detailed protocol

To access possible functional relationship between the COPS7B subunit of the COP9 signalosome and the circadian clock, 100 mm semi-confluent plates of circadian model cell line U2OS:Bmal1-luc (Maier et al., 2009) were transfected with a set of three Stealth short interfering RNAs targeting COPS7B and COPS4 (25 nM, ThermoFisher Scientific, Waltham, MA). Lipofectamine RNAiMAX transfection reagent was used for efficient siRNAs delivery and the protocol was performed according to the manufacturer’s directions (ThermoFisher Scientific). To determine the interaction between COPS7B and clock proteins, 100 mm plates of HEK293T were transfected with FLAG epitope-tagged clock proteins BMAL1, PER2 and CRY1 (5 ug) in the presence or absence of siRNAs targeting COPS7B (25 nM). The protocol for Lipofectamine 2000 Transfection Reagent (Invitrogen)-mediated transient transfection was used according to the manufacturer’s instructions.

Protein isolation, Immunoprecipitation, Western Blots and Antibodies

Preparation of nuclear protein extract

Request a detailed protocol

72 hr post-transfection, HEK293T or U2OS cells transiently transfected with FLAG epitope-tagged clock proteins (BMAL1, PER2 and CRY1) were synchronized by dexamethasone as described above, and harvested at the maximal protein levels for BMAL1 (12 hr after dexamethasone shock), PER2 and CRY1 (2 hr after dexamethasone shock) to obtain nuclear protein lysates. Nuclear extracts were prepared according to the method of Schreiber et al., 1989.

Preparation of whole cell protein extract

Request a detailed protocol

RIPA lysis buffer (0.02M Tris/pH7.2, 0.15M NaCl, 1% Trion X-100, 1% Na-deoxycholate, 0.1% SDS) was used to prepare whole cell protein extracts. Briefly, transfected HEK293T resp. U2OS cells were washed twice with ice-cold 1xPBS. 300 μl of RIPA lysis buffer supplemented with protease inhibitor cocktail and phosphatase inhibitors (1 mM NAF, 0.1 mM Na3VO4) was added to the 100 mm plates. After 10 min incubation on ice, cells were scraped vigorously and the suspension was spun down at 10000 RPM at 4°C for 10 min. The supernatant, containing proteins, was carefully transferred into the fresh microfuge tube. The concentration of the proteins was determined using the Pierce BCA Protein assay kit (ThermoFisher Scientific).

Immunoprecipitation 

Request a detailed protocol

was performed using standard procedures with the below mentioned adjustments (Ausubel, 2003). Nuclear extracts (1000 μg) were bound with antibodies (Monoclonal ANTI-FLAG M2 Antibody/Sigma, anti COPS7B antibody/Abcam, Cambridge, UK) for 90 min at 4°C, on a horizontal rotator at 10 RPM. The antibody-protein complex was then incubated for 1 hr with protein A agarose beads (Diagenode, Denville, NJ) at 4°C, on a horizontal rotator at 10 RPM. After incubation, the beads were washed once with high-salt buffer (50 mM Tris, 500 mM NaCl, 1% NP-40) followed by two gentle washes in low-salt buffer (50 mM Tris, 120 mM NaCl, 0.5% NP-40). Elution of the proteins from the beads was performed for 3 min at 95°C with 15 μl of 2x SDS sample buffer, containing β-mercaptoethanol. Equal amounts of IP reaction mixtures were loaded on a 9% SDS-PAGE gel together with 1/20 of the IP amounts of input. The protein gel electrophoresis and Western blotting was performed using standard procedures (Ausubel, 2003). Equal loading and size detection using a protein ladder were verified by Ponceau-S staining of membranes prior to probing.

Antibodies 

Request a detailed protocol

For co-immunoprecipitation, Monoclonal ANTI-FLAG M2 Antibody from Sigma (F3165 RRID:AB_259529) and anti-COPS7b antibody from Abcam (ab124718 RRID:AB_10971678) were used at 1:50 dilution. For detection in co-immunoprecipitation experiments, primary Monoclonal ANTI-FLAG M2 Antibody from Sigma was diluted at 1:1000, primary anti-COPS7b antibody from Abcam was used at 1:1000. For detection of protein levels in an ‘’around-the-clock’’ experiment primary anti-COPS7b antibody from Abcam was used at 1:1000 and Anti-βIII Tubulin mAB from Promega (G7121 RRID:AB_430874) was used at 1:500 dilutions. The probing of the secondary antibody was done at 1:10000 for IRDye 680–goat anti-mouse IgG (926–32220 RRID:AB_621840; Licor, Lincoln, NE) and 1:10000 for IRDye 800–goat anti-rabbit IgG (926–33210 RRID: AB_10796098; Licor).

Metabolic Pulse-chase

Request a detailed protocol

The half-lives of transfected FLAG epitope-tagged circadian clock proteins (BMAL1, PER2 and CRY1) in the absence (as a control) and presence of siRNAs targeting COPS7B were measured by metabolic pulse chase protein labeling according to the protocol of Zhou et al. Briefly, 72 hr post-transfection, 100 mm plates of HEK293T cells were starved of amino acids for 1 hr and metabolically labeled for 30 min with 100 μCi of Express-35S protein-labeling mix (PerkinElmer Life Sciences, Schwerzenbach, CH) in 1 ml methionine- and cysteine-free DMEM medium supplemented with 10% dialysed FCS (‘the pulse’). Sunsequently, cells were washed to remove unbound radioactive amino acids, and incubated in complete medium (‘the chase’). At indicated time points (0, 2, 4, and 6 hr), the cells were disrupted in ice-cold RIPA lysis buffer (0.02M Tris/pH7.2, 0.15M NaCl, 1% Trion X-100, 1% Na-deoxycholate, 0.1% SDS), and 1000 μg of whole cell extract were immunoprecipitated with Monoclonal ANTI-FLAG M2 Antibody from Sigma (F3165). Proteins were subjected to 9% SDS-PAGE and the gels were dried for 2 hr at 80°C by Model 583 gel dryer (BIO-RAD, Hercules, CA). Labeled proteins were visualized by autoradiography and quantitative analysis of three independent experiments for each circadian clock protein was performed with ImageStudio software.

Preparation nuclear protein extracts for SILAC MS analysis

Livers were homogenized in sucrose homogenization buffer containing 2.2 M sucrose, 15 mM KCl, 2 mM EDTA, 10 mM HEPES (pH 7.6), 0.15 mM spermin, 0.5 mM spermidin, 1 mM DTT, and protease inhibitors (0.5 mM PMSF, 10 µg/ml Aprotinin, 0.7 µg/ml Pepstatin A, and 0.7 µg/ml Leupeptin). Lysates were deposited on a sucrose cushion containing 2.05 M sucrose, 10% glycerol, 15 mM KCl, 2 mM EDTA, 10 mM HEPES (pH 7.6), 0.15 mM spermin, 0.5 mM spermidin, 1 mM DTT, and protease inhibitors. After 45 min of centrifugation at 105,000 g at 4°C. The nuclei were suspended in a nucleus buffer composed of 10 mM HEPES (pH 7.6), 100 mM KCl, 0.1 mM EDTA, 10% Glycerol, 0.15 mM spermine, 0.5 mM spermidine, 0.1 mM NaF, 0.1 mM sodium orthovanadate, 0.1 mM ZnSO4, 1 mM DTT, and protease inhibitors (0.5 mM PMSF, 10 μg/ml Aprotinin, 0.7 μg/ml Pepstatin A, and 0.7 μg/ml Leupeptin). Then nuclear protein extracts were obtained by adding an equal volume of NUN buffer (2 M urea, 600 mM NaCl, 50 mM HEPES (pH 7.6), 1 mM DTT, protease inhibitors (cOmplete ULTRA, Roche), a phosphatase inhibitor cocktail (PhosphoSTOP, Roche), and deacetylase inhibitors AGK7, salermide, and trichostatin A (all from SantaCruz Biochemicals, Dallas, TX), followed by a 20 min incubation on ice. The supernatants resulting from a 10 min centrifugation at 21,000 g at 4°C constituted the nuclear extracts (NE) and were used to perform proteins MS-analysis. Protein extracts were quantified using a BCA protein assay kit (Thermo Fisher Scientific).

Nuclear protein SILAC-MS and data analysis

Request a detailed protocol

Mixed samples were reduced with 5 mM DTT and alkylated with 18.75 mM iodoacetamide for 30 min at room temperature in the dark. After an ethanol-acetate precipitation, they were resuspended in 250 mM triethylammonium bicarbonate pH 8.0 containing 4M urea and digested overnight at 37°C with sequencing grade modified trypsin (Promega, Madison, WI) at a 1:50 (w/w) trypsin:protein ratio. The obtained peptide mixtures (250 µg total material) were desalted on SepPak C18 cartridges (Waters Corp., Milford, MA), dried, dissolved in 4M Urea with 0.1% Ampholytes pH 3–10 (GE Healthcare) and fractionated by off-gel focusing as described (Geiser et al., 2011). The 10 fractions obtained were desalted on a microC18 96-well plate (Waters Corp.), dried and resuspended in 0.1% formic acid, 3% (v/v) acetonitrile for LC-MS/MS analysis. Samples were analyzed on a hybrid linear trap LTQ-Orbitrap Velos mass spectrometer (Thermo Fisher Scientific) interfaced via a nanospray source to a Dionex RSLC 3000 nanoHPLC system (Dionex, Sunnyvale, CA). Peptides were separated on a reversed-phase Acclaim Pepmap nanocolumn (75 μm ID x 25 cm, 2.0 μm, 100 Å, (Dionex)) with a gradient from 5% to 85% acetonitrile in 0.1% formic acid (total time: 200 min) and a flow rate of 300 nl/min. Full MS survey scans were performed at 60’000 resolution. In data-dependent acquisition controlled by Xcalibur 2.1 software (Thermo Fisher Scientific), the twenty most intense multiply charged precursor ions detected in the full MS survey scan were selected for CID fragmentation in the LTQ linear trap with an isolation window of 4.0 m/z and then dynamically excluded from further selection during 35 s.

MS data were analyzed and quantified with MaxQuant version 1.3.0.5 (Cox and Mann, 2008), using Andromeda as search software (Cox et al., 2011) against UniProt (release 2012_02) database restricted to mouse (Mus musculus) taxonomy and a custom database containing usual contaminants (digestion enzymes, keratins, etc). Cleavage specificity was trypsin/P (cleavage after K, R, including KP and RP) with two missed cleavages. Mass tolerances were of 6 ppm for the precursor and 0.5 Da for CID tandem mass spectra. The iodoacetamide derivative of cysteine was specified as a fixed modification, and oxidation of methionine and protein N-terminal acetylation were specified as variable modifications. Protein identifications were filtered at 1% FDR established by MaxQuant against a reversed sequence database. A minimum of one unique peptide was necessary to discriminate sequences which shared peptides. Details of peak quantitation and protein ratio computation by MaxQuant are described elsewhere (Cox and Mann, 2008). Raw mass spectrometry data and search engine outputs have been deposited to the ProteomeXchange Consortium (proteomexchange.org) via the PRIDE partner repository (Vizcaíno et al., 2016) with the dataset identifier PXD003818.

Rhythmicity analysis for nuclear proteins

Request a detailed protocol

We assessed the rhythmicity in temporal nuclear accumulation of proteins using harmonic regression, as previously (Mauvoisin et al., 2014). Briefly, focusing on only 24 hr periodicity (data are generated under 24 hr LD cycles), we used a multiple linear regression for each relative protein time trace y(t) (transformed to log2 units). For this analysis, we used the relation:

y(t)=μ+ asin(2πt24)+ b cos(2πt24)+noise

where μ is the mean, whereas aand bare the coefficients of cosine and sine functions with period of 24 hr, respectively. The resulting p-values of all proteins are used to estimate false discovery rate (FDR) by the Benjamini–Hochberg method (Benjamini and Hochberg, 1995).

Data availability

The following previously published data sets were used
    1. Dermitzakis E
    (2011) Gencord
    Publicly available at the NCBI Sequence Read Archive (accession no. EGAD00000000027).

References

  1. Book
    1. Ausubel FM
    (editors) (2003)
    Current Protocols in Molecular Biology
    Brooklyn, New York: John Wiley.
    1. Benjamini Y
    2. Hochberg Y
    (1995)
    Controlling the false discovery rate: a practical and powerful approach to multiple testing
    Journal of the Royal Statistical Society. Series B 57:289–300.
    1. Brown SA
    2. Azzi A
    (2013)
    Handbook of Experimental Pharmacology
    66–45, Peripheral circadian oscillators in mammals, Handbook of Experimental Pharmacology.
    1. Cheong JK
    2. Virshup DM
    (2011) Casein kinase 1: Complexity in the family
    The International Journal of Biochemistry & Cell Biology 43:465–469.
    https://doi.org/10.1016/j.biocel.2010.12.004
    1. Nica AC
    2. Dermitzakis ET
    (2013) Expression quantitative trait loci: present and future
    Philosophical Transactions of the Royal Society B: Biological Sciences 368:20120362.
    https://doi.org/10.1098/rstb.2012.0362

Decision letter

  1. Patrick Nolan
    Reviewing Editor; Medical Research Council Harwell, United Kingdom

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "The genomic landscape of human cellular circadian variation points to a novel role for the human signalosome" for consideration by eLife. Your article has been reviewed by three peer reviewers, one of whom, Patrick Nolan (Reviewer #1), served as Guest Reviewing Editor, and the evaluation has been overseen by Mark McCarthy as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Yoshitaka Fukada (Reviewer #2) and Michael McCarthy (Reviewer #3).

The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.

Summary:

Gaspar et al. present their findings on circadian variation in human umbilical cord fibroblasts. They used a GWAS study and identified SNPs associated with a number of molecular circadian phenotypes. The authors claim the genetic variants they identified show greater than expected differences in allele frequency in a group of subjects with extreme chronotypes. Based on GWAS data, they focused on COPS7B, a subunit of the COP9 signalosome, for a number of mechanistic investigations into period determination. The authors use a number of approaches, including RNAi, Co-IP of FLAG-tagged proteins and SILAC, to validate their initial findings and propose some interesting and novel molecular mechanisms. Specifically, they argue that COPS7B may promote the stability of the circadian transcriptional regulator, BMAL1. Based on these observations, they posit that clock protein half-life is determined by contributions both from destabilising factors via the ubiquitin proteasome and stabilising factors via the signalosome. The latter would represent a significant addition to our understanding of the molecular basis of circadian period determination.

Essential revisions:

GWAS:

Conducting a GWAS study with such a small number of samples (140-147) is ambitious and this probably underlies the high p-values for the data described throughout the manuscript. Although such datasets usually require higher power, the reviewers understand the limitations here. It is not reasonable to expect numbers in a cell-based experiment to go much larger, or equal to those seen in GWAS (N>10,000). Both points should be discussed and justified.

Please justify multiple testing using a threshold of FDR 0.1. Why choose this threshold? An FDR of q<0.05 is reasonable for a GWAS study, even a cell based study. An FDR of 0.1 seems unusually high. If not justifiable, then perhaps gene list should be modified.

In subsequent text, the authors quote the p value for the two most significant SNPs affecting PER2 expression, not the q value. Why?

It would be helpful to know the degree of phenotypic correlation across the 5 rhythm parameters, and to estimate their genetic overlap. Do any of the markers associate with more than one parameter (as expected)? Or if not, why? For example, in comparing Manhattan plots, PER2 expression is the only variable where the Chr2 SNPs show anything near a threshold. Given the eventual mechanistic validation (BMAL1 stability), shouldn't we expect a similar pattern in the Manhattan plot for PER1 expression? Please provide this information and justify.

In Figure 1, the statistical method is not well described, nor is a citation provided that validates the approach. The reviewers have concerns that simply calculating the difference in allele frequency between groups may not be valid across the entire range of allele frequencies. Why wasn't a more standard approach employed (e.g. chi square analysis, or for selected SNPs, trait morningness vs. genotype association)? If possible, these analyses should be included. Also in Figure 1, is it appropriate to study all association with chronotype equally? For instance, SNP associations with phase may be predicted to have stronger contribution to chronotype than amplitude. This approach should be justified, especially in light of the small sample size and the multiple comparisons problem.

RNAi:

Aspects of the RNAi experiments are hard to interpret. Of the 28 genes affecting period and 29 affecting amplitude, is there overlap (i.e. genes that affect both)? Of the control RNAis, the results should more clearly state how many control RNAis were used, and what proportion affected period or amplitude. Does a chi-square test indicate this proportion is significantly different?

In subsection “Genes associated with positive SNPs affect circadian clock function”, when assessing the probability that the genes affect the circadian clock when compared to random genesets, it is unclear why the authors have used one test to assess the reduced & increased period/amplitude and a different test for the combined cumulative distribution. Presumably this is from analysis of the same data so the same tests should be used, please change or justify.

Gene enrichment analysis, Figure 3:

The authors state that they 'applied gene ontology enrichment analysis to the list that we obtained'. There is a bit of confusion here. Was this applied to the 59 genes from 3C or the 76 genes from the original eQTL? What threshold was used for the GSEA? Figure 3A and 3B should be amalgamated into a graph with clear p-values added to the graph. As written in subsection “Positive SNPs highlight protein catabolism as essential to human circadian variation”, paragraph 2, authors then go on to pick the genes from 5 of the top 6 pathways for further analysis shown in Figure 3D-F. Is this somehow relevant to what is shown in Figure 3C? Only some of the pathways shown in Figure 3C are significant and included in the top 6 pathways. For example, oxidoreductase activity (included) is not in the top 6 pathways in Figure 3—figure supplement 1 while macromolecule catabolic process (not included) is? Finally, given that the resulting genes are not significant in the GWAS and are not from a random sample, please provide evidence that the significant increase/decrease in period (Figure 3E) is not by chance.

Figure 4 and associated text:

Is the rs920400 AA genotype associated with long or short period in fibroblasts? How can this be reconciled with 'average PER2 expression' being higher? Is it justifiable to use a single timepoint here to look at PER2 expression? Do the genetic association findings of COPS7B share any association with the BMAL-Luc parameters? As the main result, the COPS7B association should be reported for all of the 5 parameters selected. If the genetic association does not match the RNAi experiment, the discrepancy should be discussed. It is difficult to accept that COPS4 knockdown lengthens period (by 9.9 hrs!) as it appears to completely obliterate the bioluminescence rhythm. What is the condition of the cells after knockdown? Labelling in Figure (COP4) is incorrect.

Figure 5 and associated text:

The authors claim direct binding of COP9 signalosome with several clock proteins, but the data need to be polished in order to fully support this conclusion. The immunoprecipitation analysis (Figure 5A, B, D) lacks essential experimental controls: in the FLAG-IP experiment, COPS7b signal should be no longer be detected in the absence of co-expression of FLAG-proteins in the same gel. These controls should be included in Western Blots. In the FLAG IP (a) the COPS7B antibody appears to cross-react with both PER2-FLAG and CRY1-FLAG. This confound should be addressed and justified. Both panels seem to show that interaction of COPS7B with PER2-FLAG and CRY1-FLAG is greater than with BMAL1-FLAG, yet the authors are proposing a direct interaction with BMAL1, how is this so? Authors should address this as it affects their arguments in the paper. Figure 5C. What does the * signify and how was it measured? Presumably this should be a repeated measures ANOVA determination. Statistical methods should be included. How does this translate into 'the half-life of BMAL1 was reduced twofold'? Figure 5—figure supplement 1. Shouldn't you be seeing some protein degradation in these graphs, there is nothing evident here, particularly so in the CRY1 experiment. Some explanation is required in the text. How is this related to 'The half lives of CRY1 and PER2 remained unchanged'? Neither levels approach anything near 50% of the time-0 timepoint. What does 'even if their absolute abundance varied' mean?

Subsection “COPS7B and the COP9 signalosome are rhythmically imported into the nucleus”:

The argument the authors make about the conflicting phases of expression of PER/CRY with COPS7B in U2OS cells might hold for native protein interactions but this should not hold when constitutively ovexpressing recombinant PER-FLAG and CRY-FLAG proteins. There must be some other cause for the discrepancy seen in protein interactions in the two cell lines. Please address these issues in Results and Discussion.

Presentation and style:

There are numerous instances of errors in cross-referencing of Figures, referencing citations and overall flow of text in Results section.

Abstract:

Given that the Abstract is the first port-of-call for many readers, it should be a realistic reflection of the main findings and conclusions in the manuscript. Some text such as 'Overwhelmingly gene set enrichment points to differences in protein catabolism as the major source of clock variation in humans' is rather strong. The reviewers request that this be toned down appropriately.

https://doi.org/10.7554/eLife.24994.027

Author response

Essential revisions:

GWAS:

Conducting a GWAS study with such a small number of samples (140-147) is ambitious and this probably underlies the high p-values for the data described throughout the manuscript. Although such datasets usually require higher power, the reviewers understand the limitations here. It is not reasonable to expect numbers in a cell-based experiment to go much larger, or equal to those seen in GWAS (N>10,000). Both points should be discussed and justified.

These are indeed a valid point, and we thank the reviewers for their understanding of the problems of cell-based genetics. (Even as it stands, ours was a story close to a decade in the making.) We have now discussed the points raised by this reviewer:

“Typically, because of the labor involved in collecting and analyzing cellular material, the cohorts used in such investigations are much smaller (hundreds of samples rather than >10,000 in conventional GWAS). Partially compensating for this deficiency is the relative simplicity of the analytical system and the precision of expression trait measurement.”

and again:

“While multiple suggestive alleles were identified, it should be noted that none achieved “genome wide” significance, as might be expected from a cell-based study of moderate sample size and 2.5million SNPs tested. Therefore, we undertook further experiments to demonstrate the relevance of our results.”

Please justify multiple testing using a threshold of FDR 0.1. Why choose this threshold? An FDR of q<0.05 is reasonable for a GWAS study, even a cell based study. An FDR of 0.1 seems unusually high. If not justifiable, then perhaps gene list should be modified.

Given the multistep selection process we used, we preferred to tolerate a larger rate of false positives initially (1/10 vs. 1/20) that would be screened out later. It should be noted that “suggestive” significance for GWAS is generally defined at 10e-5, and this corresponds in our study to an FDR of 0.1. It should further be noted that our FDR is based upon a strict correction for the 2.5million SNPs tested, and for European genomes the 2.5M SNP chip that we used is in fact “far too dense”: it is estimated that only 260,000 “tag-SNPs” are sufficient to capture variance in all common SNPs in the Phase I Hapmap project for European genomes. Thus, a Bonferroni-based correction is in reality 10x too severe! An excellent discussion of this point can be found here: Nature. 2005 Oct 27; 437(7063): 1299– 1320.

Finally, the validation of this strategy is in our opinion visible by the success of the GSEA approach, which demonstrated the functionally related significance of these hits.

We now discuss our rationale explicitly in the text:

“We chose this threshold to correspond to “suggestive” significance (uncorrected p<1x10-5) in genome-wide association: since it is estimated that only a tenth of the SNPs on the high-density chip that we used would be necessary to capture all common variation in individuals of European origin, a corrected FDR is therefore up to tenfold “over-corrected”. We therefore preferred to select genes at a relaxed FDR, and subsequently apply further criteria.”

In subsequent text, the authors quote the p value for the two most significant SNPs affecting PER2 expression, not the q value. Why?

For the “raw” GWAS, we wanted the values to be immediately comparable across other genetic studies, and p values have been reported much longer than q values. The main reason for the persistence in the use of p-values in GWAS is that correction for multiple testing after genotyping with high-density arrays is highly disputable. Essentially, the 2.5million SNPs on a modern chip are not independent. In fact, the Hapmap project shows that the oversampling rate is tenfold for Europeans. By leaving these values as p-values, anyone can easily correct however they want.

It would be helpful to know the degree of phenotypic correlation across the 5 rhythm parameters, and to estimate their genetic overlap. Do any of the markers associate with more than one parameter (as expected)? Or if not, why? For example, in comparing Manhattan plots, PER2 expression is the only variable where the Chr2 SNPs show anything near a threshold. Given the eventual mechanistic validation (BMAL1 stability), shouldn't we expect a similar pattern in the Manhattan plot for PER1 expression? Please provide this information and justify.

We now include the Spearman correlation pairwise across all phenotypes as a way of estimating genetic overlap. We observe varying degrees of global correlation across parameters – i.e. a nonrandom association of p values across the different phenotypes assayed. As might be expected, in particular PER2 expression correlates with PER1 expression, and period correlates with phase (a new Table 1). However, this reviewer correctly notes that correlation across top hits is surprisingly absent – i.e. wouldn’t a top hit for PER2 expression also be one for PER1 or period length? Apparently not.

We now note this explicitly in the text:

“We next compared SNPs across the traits that we identified. Considering only alleles achieving “suggestive” genome-wide significance (p<10-5), we see correlation between multiple pairs of traits, notably PER1 and PER2 expression, and period and phase (Table 1). Given the close homology and overlapping function of the PER proteins, and the known correlation between period length and phase, these correlations are expected. Surprisingly, however, the most significant alleles in any one category are not among the most significant alleles in another, suggesting relatively independent genetic regulation of key circadian clock parameters (“state variables”) for small changes, and the resilience of the overall circadian mechanism to small perturbations.”

We continue the theme in our Discussion:

“It is worth noting that although knockdown of COP9 subunits altered circadian period, phase, and PER2 expression, the COPS7B-associated alleles rs920400 and rs10195385 that we genetically identified were associated significantly only with PER2 expression in vitro. More generally, although there was detectable and significant genotypic correlation between period and phase, and PER1 and PER2 expression in our screen, global genetic correlation among these different phenotypes was surprisingly weak. We attribute these differences to the idea that small changes in cellular clock properties can be robustly tolerated without affecting other clock parameters. Such an idea has also been supported by cellular dosage studies individual clock proteins.”

In Figure 1, the statistical method is not well described, nor is a citation provided that validates the approach.

We thought long and hard about how to demonstrate an enrichment that would be valid in a moderately-sized pool of only 40 individuals of extreme chronotype, and invented this method for the task. The method is analogous to those that verify individual eQTL hits in biochemical pathways by validating unequal allele penetrance in individuals with diseases affecting that pathway (e.g. Loeuillet et al., 2008).

We now describe it in better detail in a new section of the Materials and methods.

The reviewers have concerns that simply calculating the difference in allele frequency between groups may not be valid across the entire range of allele frequencies.

This is an important issue that we considered carefully – e.g. an extremely rare allele, present in only one subject who was an extreme early or late type, would skew the results at stringent p values where relatively few alleles were present. We combated this in multiple ways: first, our study already removed rare alleles (minor allele <5%) from consideration. Secondly, we included these data in two different versions: one with all SNPs, and one with common tag-SNPs (Figure 1B, C and Figure 1—figure supplement 4). This corrects for the multiple comparison problem – the tag-SNP approach counts all alleles in a hapblock only once – and also by definition eliminates very rare alleles). Finally, to fully satisfy this concern, we have gone back and looked specifically at all SNPs p<10-6. The rarest minor allele among all these SNPs is 0.1524. Therefore, we believe that this concern, though valid, is not a problem in our analysis.

Why wasn't a more standard approach employed (e.g. chi square analysis, or for selected SNPs, trait morningness vs. genotype association)? If possible, these analyses should be included.

We respectfully disagree with this idea. Because of the limited number of extreme chronotypes in our pool, we do not think that it is fair to pick single alleles and show their enrichment. This is something that would be easy to do, but we do not think it would be valid. The whole strength of the method we devised is that it avoids singling out any particular allele, and shows that on the whole alleles considered at more stringent p values have greater influence upon chronotype than those at lesser.

Also in Figure 1, is it appropriate to study all association with chronotype equally? For instance, SNP associations with phase may be predicted to have stronger contribution to chronotype than amplitude. This approach should be justified, especially in light of the small sample size and the multiple comparisons problem.

Given the limited number of alleles we obtain, we think that it is important to consider them collectively in order to add power to subsequent pathway analyses. Moreover, while it is commonly predicted that some clock properties affect chronotype more than others, this relationship is far from clear in the literature. For example, we ourselves have shown that amplitude directly affects entrained phase, because a clock of high amplitude will show a smaller phase shift response. Thus, humans with high amplitude tend to be later chronotypes (Brown, PNAS 2008).

RNAi:

Aspects of the RNAi experiments are hard to interpret. Of the 28 genes affecting period and 29 affecting amplitude, is there overlap (i.e. genes that affect both)?

At the biochemical level, absolutely – COP9 subunits are one particularly relevant example. At the genetic level, mostly not – i.e., as already mentioned in the context of Figure 1, we do not find highly significant polymorphisms affecting both parameters. However, we do not find this so surprising: in vitro, black-and- white changes in one circadian parameter definitely affect another. However, with smaller changes, these parameters behave much more independently. A good example of this can be seen in the dose-dependent RNAi against clock genes seen in Baggs et al., 2009.

Of the control RNAis, the results should more clearly state how many control RNAis were used, and what proportion affected period or amplitude. Does a chi-square test indicate this proportion is significantly different?

We believe that a misunderstanding occurred here, and we have therefore rewritten our description of this experiment in subsection “Positive SNPs highlight protein catabolism as essential to human circadian variation” and in the corresponding figure legends. Although only a subset of randomly chosen control RNAis were shown as discrete black circles, the underlying curves represented the entire genome, and the statistical tests described were against this control. However, it is impossible to plot aesthetically so many individual values. Statistical tests for all comparisons except one are significant (as noted in the text), and all exact p-values are cited in the text.

In subsection “Genes associated with positive SNPs affect circadian clock function”, when assessing the probability that the genes affect the circadian clock when compared to random genesets, it is unclear why the authors have used one test to assess the reduced & increased period/amplitude and a different test for the combined cumulative distribution. Presumably this is from analysis of the same data so the same tests should be used, please change or justify.

We have now homogenized our presentation here (Figures 2B, C and 3E, F). In fact, both types of analyses were already present in both comparisons (against the whole genome, not random genesets as mentioned above), but one was in the main figure and the other in supplements. We chose two different representations so that the reader could see in one case the change in magnitude of phenotype, and in the other case a directional shift in the distribution. That is also why two different statistical tests were cited. We think that this approach is valid, but in any case both approaches and both tests were used in both cases, and now the figures look homogenous.

Gene enrichment analysis, Figure 3:

The authors state that they 'applied gene ontology enrichment analysis to the list that we obtained'. There is a bit of confusion here. Was this applied to the 59 genes from 3C or the 76 genes from the original eQTL?

Figures 3A and 3B were generated from the original list obtained purely by individual gene q-value. 3c-f were a completely separate analysis done using the GSEA algorithm, and we highlight the similarity of the results in the text. Thus, these are two completely separate analyses using different methods that nevertheless obtained similar results.

We have added additional phrases to make this distinction clearer.

“We first applied gene ontology enrichment analysis to the list of 60 genes…”

What threshold was used for the GSEA?

GSEA intentionally uses a highly relaxed stringency for individual alleles (10e-4) and then specifically looks for enriched sets. The thresholds we have used here are consistent with those used previously, and recommended by the authors of the tool. We have now noted this relaxed stringency explicitly in the subsection “COPS7B influences PER2 expression and clock function”.

Figure 3A and 3B should be amalgamated into a graph with clear p-values added to the graph.

We have now exchanged the pathway diagram of Figure 3A for a simplified graph, and added as supplements the full pathway diagram, as well as a supplementary file with IDs for the GO biological processes, the genes from our gene set within the pathways, and the adjusted p-values. Figure 3B, the GOslim functional distribution, represents complementary information. It shows the highest-level functional annotations of the groups of genes identified. Since it is showing the actual annotation rather than the probability of mapping by chance, a p value is not appropriate, and indeed is not supplied with the tool.

As written in subsection “Positive SNPs highlight protein catabolism as essential to human circadian variation”, paragraph 2, authors then go on to pick the genes from 5 of the top 6 pathways for further analysis shown in Figure 3D-F. Is this somehow relevant to what is shown in Figure 3C? Only some of the pathways shown in Figure 3C are significant and included in the top 6 pathways. For example, oxidoreductase activity (included) is not in the top 6 pathways in Figure 3—figure supplement 1 while macromolecule catabolic process (not included) is?

Clearly, this was not very well explained. The raw GSEA output is shown in Figure 3—figure supplement 2. Five of the top six pathways were different aspects of protein catabolism, and this is what we chose to investigate in the rest of the paper. Thus, the following panels d-f, and the remainder of the paper, deal only with protein stability. However, several other pathways were identified at p<0.003 or better, and we listed each only once in the main Figure 3. We now explain this abbreviation better in the text:

“The full list from GSEA analysis is shown in Figure 3—figure supplement 2A. We also show an abbreviated list, eliminating major redundancies in pathway annotations, in Figure 3C. Since by far the most significant category identified was protein catabolism, we chose to explore further the significance of alleles in this group.”

Finally, given that the resulting genes are not significant in the GWAS and are not from a random sample, please provide evidence that the significant increase/decrease in period (Figure 3E) is not by chance.

For this, we used a Kolmogorov-Smirnov test to show that the distributions are different. The comparison is against the genome as a whole, and the KS test p- value already accounts for differences in the number of elements in the two distributions.

Figure 4 and associated text:

Is the rs920400 AA genotype associated with long or short period in fibroblasts? How can this be reconciled with 'average PER2 expression' being higher?

Short period. There is an obvious contradiction here that we already discuss in the first paragraph of the Discussion. Deletion studies by multiple other labs show that PER protein knockdown is associated with a short period. We think that the increase in Per2 RNA levels is in fact an indirect effect due to an increase in BMAL1 protein.

Is it justifiable to use a single timepoint here to look at PER2 expression?

While many timepoints might be useful for a circadian question, it was not feasible within the context of the hundreds of samples we examined. Since we attribute the transcriptional phenotype to an indirect effect, we decided not to pursue Per2 expression in further detail. Clearly, altered timing of expression would be reflected in changed levels at a single timepoint. However, our main goal here was to verify the findings of the screen itself.

We have added the following text so that the reader is aware of this potential confound:

“It should be noted that here and in the screen itself, a single timepoint was used to investigate PER2 expression, corresponding to its expected first peak after cell line synchronization. Thus for any single subject, absolute changes in level could also reflect changes in timing. Therefore, we next explored this effect using exogenous assays.”

Do the genetic association findings of COPS7B share any association with the BMAL-Luc parameters? As the main result, the COPS7B association should be reported for all of the 5 parameters selected. If the genetic association does not match the RNAi experiment, the discrepancy should be discussed.

The only striking association is with PER2 expression. We now report all associations as requested. Like this reviewer, I would have preferred that a striking association with PER2 expression would be associated with an equally striking change in period. This is not the case. As requested, we now discuss this explicitly in the text:

“It is worth noting that although knockdown of COP9 subunits altered circadian period, phase, and PER2 expression, the COPS7B-associated alleles rs920400 and rs10195385 that we genetically identified were associated significantly only with PER2 expression in vitro. More generally, although there was detectable and significant genotypic correlation between period and phase, and PER1 and PER2 expression in our screen, global genetic correlation among these different phenotypes was surprisingly weak. We attribute these differences to the idea that small changes in cellular clock properties can be robustly tolerated without affecting other clock parameters. Such an idea has also been supported by cellular dosage studies individual clock proteins.”

It is difficult to accept that COPS4 knockdown lengthens period (by 9.9 hrs!) as it appears to completely obliterate the bioluminescence rhythm. What is the condition of the cells after knockdown? Labelling in Figure (COP4) is incorrect.

This reviewer refers to the fact that cellular toxicity lengthens circadian period. However, these cells were viable and continued to grow, but are basically arrhythmic. The reduction in bioluminescence was not due to cell death. We have highlighted this effect better by now including raw bioluminescence traces.

Figure 5 and associated text:

The authors claim direct binding of COP9 signalosome with several clock proteins, but the data need to be polished in order to fully support this conclusion. The immunoprecipitation analysis (Figure 5A, B, D) lacks essential experimental controls: in the FLAG-IP experiment, COPS7b signal should be no longer be detected in the absence of co-expression of FLAG-proteins in the same gel. These controls should be included in Western Blots.

Since we did the immunoprecipitation in both directions, we initially judged this control redundant. However, we have now done this control as suggested by the reviewer and we see no false positive signal. The results from this experiment are now included in Figure 5 —figure supplement 1.

In the FLAG IP (a) the COPS7B antibody appears to cross-react with both PER2-FLAG and CRY1-FLAG. This confound should be addressed and justified.

This is not the case, but is rather an unfortunate background of slightly different size. Short of trimming the bands more closely, we can’t make this problem “disappear”. We now make reference to it explicitly in the legend.

Both panels seem to show that interaction of COPS7B with PER2-FLAG and CRY1-FLAG is greater than with BMAL1-FLAG, yet the authors are proposing a direct interaction with BMAL1, how is this so? Authors should address this as it affects their arguments in the paper.

This is an important point that we address in Discussion paragraph three, and have now expanded here and in paragraph four. In transfection experiments into noncircadian cells, we do indeed see interaction with all three-clock proteins mentioned. In circadian cells, we see only the BMAL1 interaction, and subsequently show that the COPS9 signalosome is imported into the nucleus roughly coincident with BMAL1 in cells and tissues. Thus, whether or not it could interact with all three, functionally it appears to interact primarily with BMAL1.

Figure 5C. What does the * signify and how was it measured? Presumably this should be a repeated measures ANOVA determination. Statistical methods should be included.

We were simpler here: this was only a Student T test for the replicates at each timepoint, p=0.049. We now describe it in the legend.

How does this translate into 'the half-life of BMAL1 was reduced twofold'? Figure 5—figure supplement 1. Shouldn't you be seeing some protein degradation in these graphs, there is nothing evident here, particularly so in the CRY1 experiment. Some explanation is required in the text. How is this related to 'The half lives of CRY1 and PER2 remained unchanged'? Neither levels approach anything near 50% of the time-0 timepoint. What does 'even if their absolute abundance varied' mean?

We thank the reviewer for pointing out this problem with our quantification. As the blots themselves show, degradation is indeed occurring, and it is not different between wildtype and knockdown. Values approach 50% at 2-4 hours. Curiously, though, they subsequently “flatten” during the 6 hours analysed; these experiments were repeated multiple times. While we do not understand why further degradation does not occur – this might even be an interesting regulatory phenomenon – we prefer not to speculate, and simply use the assay to calculate differences between wildtype and knockdown.

Subsection “COPS7B and the COP9 signalosome are rhythmically imported into the nucleus”:

The argument the authors make about the conflicting phases of expression of PER/CRY with COPS7B in U2OS cells might hold for native protein interactions but this should not hold when constitutively ovexpressing recombinant PER-FLAG and CRY-FLAG proteins. There must be some other cause for the discrepancy seen in protein interactions in the two cell lines. Please address these issues in Results and Discussion.

We have now qualified our conclusions as suggested by this reviewer:

“It should be noted that this result does not completely explain why interactions with other clock proteins would not be observed when they are exogenously transfected, and suggests therefore that time-specific cofactors might be required for the interaction.”

…and then:

“Although we could detect interactions with multiple transfected clock proteins in nonrhythmic HEK cells, in rhythmic U2OS cells only BMAL1 interaction was detectable, suggesting further that other time-of-day-specific factors might be required for this interaction.”

Presentation and style:

There are numerous instances of errors in cross-referencing of Figures, referencing citations and overall flow of text in Results section.

We apologise for these errors. Most result from a reformatting for eLife, followed by accidental download of the wrong version. Still, we should have noticed. We believe that all these errors have now been corrected.

Abstract:

Given that the Abstract is the first port-of-call for many readers, it should be a realistic reflection of the main findings and conclusions in the manuscript. Some text such as 'Overwhelmingly gene set enrichment points to differences in protein catabolism as the major source of clock variation in humans' is rather strong. The reviewers request that this be toned down appropriately.

We have weakened the wording to: “Gene set enrichment points to differences in protein catabolism as one major source of clock variation in humans.”

https://doi.org/10.7554/eLife.24994.028

Article and author information

Author details

  1. Ludmila Gaspar

    Institute of Pharmacology and Toxicology, University of Zurich, Zurich, Switzerland
    Present address
    Max-Planck Institute, Tuebingen, Germany
    Contribution
    Conceptualization, Investigation, Methodology, Writing—original draft, Project administration
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6673-7492
  2. Cedric Howald

    1. Department of Genetic Medicine and Development, University of Geneva, Geneva, Switzerland
    2. Institute of Genetics and Genomics in Geneva, University of Geneva, Geneva, Switzerland
    Contribution
    Data curation, Investigation, Writing—original draft
    Competing interests
    No competing interests declared
  3. Konstantin Popadin

    Department of Genetic Medicine and Development, University of Geneva, Geneva, Switzerland
    Present address
    Immanuel Kant Baltic Federal University, Kaliningrad, Russia
    Contribution
    Data curation, Investigation, Methodology
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-2117-6086
  4. Bert Maier

    Charité–Universitätsmedizin, Laboratory of Chronobiology, Berlin, Germany
    Contribution
    Resources, Investigation, Methodology, Writing—original draft
    Competing interests
    No competing interests declared
  5. Daniel Mauvoisin

    Department of Pharmacology and Toxicology, University of Lausanne, Lausanne, Switzerland
    Present address
    Institute of Bioengineering, School of Life Sciences, Ecole PolytechniqueFédérale de Lausanne and Swiss Institute of Bioinformatics, Lausanne, Switzerland
    Contribution
    Resources, Data curation, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0571-0741
  6. Ermanno Moriggi

    Institute of Pharmacology and Toxicology, University of Zurich, Zurich, Switzerland
    Contribution
    Conceptualization, Formal analysis, Investigation
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-4600-5777
  7. Maria Gutierrez-Arcelus

    1. Department of Genetic Medicine and Development, University of Geneva, Geneva, Switzerland
    2. Institute of Genetics and Genomics in Geneva, University of Geneva, Geneva, Switzerland
    Contribution
    Resources, Investigation, Methodology
    Competing interests
    No competing interests declared
  8. Emilie Falconnet

    1. Department of Genetic Medicine and Development, University of Geneva, Geneva, Switzerland
    2. Institute of Genetics and Genomics in Geneva, University of Geneva, Geneva, Switzerland
    Contribution
    Investigation, Methodology
    Competing interests
    No competing interests declared
  9. Christelle Borel

    1. Department of Genetic Medicine and Development, University of Geneva, Geneva, Switzerland
    2. Institute of Genetics and Genomics in Geneva, University of Geneva, Geneva, Switzerland
    Contribution
    Conceptualization, Resources, Data curation, Writing—review and editing
    Competing interests
    No competing interests declared
  10. Dieter Kunz

    Institute of Physiology, Charité-Universitätsmedizin Berlin, Working Group Sleep Research & Clinical Chronobiology, Berlin, Germany
    Contribution
    Conceptualization, Resources, Investigation, Methodology
    Competing interests
    No competing interests declared
  11. Achim Kramer

    Charité–Universitätsmedizin, Laboratory of Chronobiology, Berlin, Germany
    Contribution
    Conceptualization, Resources, Funding acquisition, Writing—review and editing
    Competing interests
    No competing interests declared
  12. Frederic Gachon

    Department of Pharmacology and Toxicology, University of Lausanne, Lausanne, Switzerland
    Present address
    Diabetes and Circadian Rhythms department, Nestlé Institute of HealthSciences, Lausanne, Switzerland
    Contribution
    Conceptualization, Resources, Funding acquisition, Methodology, Writing—review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9279-9707
  13. Emmanouil T Dermitzakis

    1. Department of Genetic Medicine and Development, University of Geneva, Geneva, Switzerland
    2. Institute of Genetics and Genomics in Geneva, University of Geneva, Geneva, Switzerland
    Contribution
    Conceptualization, Resources, Funding acquisition, Writing—review and editing
    Competing interests
    No competing interests declared
  14. Stylianos E Antonarakis

    1. Department of Genetic Medicine and Development, University of Geneva, Geneva, Switzerland
    2. Institute of Genetics and Genomics in Geneva, University of Geneva, Geneva, Switzerland
    Contribution
    Conceptualization, Resources, Funding acquisition, Writing—review and editing
    Competing interests
    No competing interests declared
  15. Steven A Brown

    Institute of Pharmacology and Toxicology, University of Zurich, Zurich, Switzerland
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Project administration, Writing—review and editing
    For correspondence
    Steven.brown@pharma.uzh.ch
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-5511-568X

Funding

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung (CRSII3_160741)

  • Steven A Brown

Zurich Hospital (CRPPSleep&Health)

  • Steven A Brown

Velux Fonden (923)

  • Steven A Brown

European Research Council (ERC-2010-StG-260988)

  • Frederic Gachon

Fondation Leenaards (Grant)

  • Frederic Gachon

Immanuel Kant Baltic University (5 Top 100 Russian Academic Excellence Project)

  • Konstantin Popadin

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

LG has been supported by the Swiss National Science Foundation. This work has received further support via SAB from the Swiss National Science Foundation, the Velux Foundation, and the Clinical Research Priority Program ‘Sleep and Health’ of the University of Zürich. LG and SAB are members of the Zurich Neurozentrum, a division of the Life Sciences Zürich graduate program. We thank Juergen Ripperger and Urs Albrecht for generous donation of antibodies. Research in the FG laboratory received funding from the European Research Council (through individual Starting Grants ERC-2010-StG-260988) and the Leenaards Foundation. KP was supported by the 5 Top 100 Russian Academic Excellence Project at the Immanuel Kant Baltic Federal University.

Ethics

Human subjects: All human samples used in these studies were obtained after approval of all protocols and procedures by the relevant responsible authorities (of the University Hospital Geneva, CH; and Charite Universitätsmedezin, Berlin, DE), and prior written informed consent was obtained from all subjects or their legal guardians.

Animal experimentation: All animal experiments were conducted with the approval of relevant cantonal veterinary authorities in Switzerland, after prior review of all procedures and planned experiments.

Reviewing Editor

  1. Patrick Nolan, Medical Research Council Harwell, United Kingdom

Publication history

  1. Received: January 10, 2017
  2. Accepted: September 1, 2017
  3. Accepted Manuscript published: September 4, 2017 (version 1)
  4. Version of Record published: September 15, 2017 (version 2)

Copyright

© 2017, Gaspar et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,174
    Page views
  • 295
    Downloads
  • 3
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Cancer Biology
    2. Cell Biology
    Linda Julian et al.
    Research Article

    Apoptosis is characterized by profound morphological changes, but their physiological purpose is unknown. To characterize the role of apoptotic cell contraction, ROCK1 was rendered caspase non-cleavable (ROCK1nc) by mutating Aspartate 1113, which revealed that ROCK1 cleavage was necessary for forceful contraction and membrane blebbing. When homozygous ROCK1nc mice were treated with the liver-selective apoptotic stimulus of diethylnitrosamine, ROCK1nc mice had more profound liver damage with greater neutrophil infiltration than wild-type mice. Inhibition of the damage associated molecular pattern protein HMGB1 or signalling by its cognate receptor TLR4 lowered neutrophil infiltration and reduced liver damage. ROCK1nc mice also developed fewer diethylnitrosamine-induced hepatocellular carcinoma (HCC) tumours, while HMGB1 inhibition increased HCC tumour numbers. Thus, ROCK1 activation and consequent cell contraction are required to limit sterile inflammation and damage amplification following tissue-scale cell death. Additionally, these findings reveal a previously unappreciated role for acute sterile inflammation as an efficient tumour suppressive mechanism.

    1. Cell Biology
    2. Stem Cells and Regenerative Medicine
    Emma Carley et al.
    Research Article Updated

    While the mechanisms by which chemical signals control cell fate have been well studied, the impact of mechanical inputs on cell fate decisions is not well understood. Here, using the well-defined system of keratinocyte differentiation in the skin, we examine whether and how direct force transmission to the nucleus regulates epidermal cell fate. Using a molecular biosensor, we find that tension on the nucleus through linker of nucleoskeleton and cytoskeleton (LINC) complexes requires integrin engagement in undifferentiated epidermal stem cells and is released during differentiation concomitant with decreased tension on A-type lamins. LINC complex ablation in mice reveals that LINC complexes are required to repress epidermal differentiation in vivo and in vitro and influence accessibility of epidermal differentiation genes, suggesting that force transduction from engaged integrins to the nucleus plays a role in maintaining keratinocyte progenitors. This work reveals a direct mechanotransduction pathway capable of relaying adhesion-specific signals to regulate cell fate.