Novel transcriptional networks regulated by CLOCK in human neurons.

The molecular mechanisms underlying human brain evolution are not fully understood; however, previous work suggested that expression of the transcription factor CLOCK in the human cortex might be relevant to human cognition and disease. In this study, we investigated this novel transcriptional role for CLOCK in human neurons by performing chromatin immunoprecipitation sequencing for endogenous CLOCK in adult neocortices and RNA sequencing following CLOCK knockdown in differentiated human neurons in vitro. These data suggested that CLOCK regulates the expression of genes involved in neuronal migration, and a functional assay showed that CLOCK knockdown increased neuronal migratory distance. Furthermore, dysregulation of CLOCK disrupts coexpressed networks of genes implicated in neuropsychiatric disorders, and the expression of these networks is driven by hub genes with human-specific patterns of expression. These data support a role for CLOCK-regulated transcriptional cascades involved in human brain evolution and function.

Molecular evolution on the human lineage is thought to have arisen primarily through modifications of gene regulation and expression rather than through changes in protein structure (King and Wilson 1975; Fontenot and Konopka 2014). In particular, the dramatic changes in brain-relevant phenotypes and disorders in humans compared with related species with conserved protein-coding genomes, such as chimpanzee, supports modifications of gene expression in human brains. Comparative brain genomics studies have identified hundreds of expression changes that are specific to human brains. While many of these changes are in genes involved in important brain functions such as synaptic plasticity and dendritic morphology (Khaitovich et al. 2004; Uddin et al. 2004; Liu et al. 2012), recent work has identified a specific increase in the transcription factor CLOCK in the human cortex (Babbitt et al. 2010; Konopka et al. 2012a). Furthermore, CLOCK was discovered to be a major “hub” gene in a human-specific frontal pole module compared with other primates using a coexpression network approach (Konopka et al. 2012a).

This result was surprising, as CLOCK has been studied primarily in regard to its role in circadian rhythms, a conserved feature of many organisms used to coordinate behavioral, physiological, and metabolic functions with a daylong cyclical environment (Takahashi 2017). These rhythms are observed across both prokaryotic and eukaryotic kingdoms, including cyanobacteria, fungi, insects, and mammals (Menaker et al. 1997; Dvornyk et al. 2003). Almost every mammalian cell genetically encodes a molecular clock capable of maintaining internal phases in the absence of external cues (Nagoshi et al. 2004; Welsh et al. 2004; Yoo et al. 2004). Thus, the potentially human-specific role of CLOCK in the human cortex required further investigation, as it seemed unlikely that this novel expression pattern would be related to a highly conserved circadian feature (Konopka et al. 2012a). In addition, few genes identified through comparative brain genomics studies have been functionally investigated. The functional study of individual genes identified from comparative genomics is therefore important to validate the increasing number of differentially expressed genes (DEGs) identified using constantly improving high-resolution next-generation sequencing technologies (Fontenot and Konopka 2014).

The transcriptional landscape downstream from CLOCK has been described previously in the mouse liver and pancreas, where CLOCK was shown through chromatin immunoprecipitation sequencing (ChIP-seq) to bind to several thousand sites as part of a heterodimeric complex with its dimerization partner, ARNTL (Yoshitane et al. 2009; Koike et al. 2012; Perelis et al. 2015). This heterodimer binds primarily at consensus E-box DNA motifs and thereafter recruits other transcription factors to impact gene expression (Ripperger and Schibler 2006). In vitro studies in 293T (Alhopuro et al. 2010) and THP-1 cells derived from an acute monocytic leukemia patient (Puram et al. 2016) have also described CLOCK binding throughout the human genome in nonneuronal systems.

CLOCK has been implicated in a number of disease states outside of the core circadian system, the most well-characterized of which are metabolic diseases. The dominant-negative mutant Clock-Δ19 mice (Vitaterna et al. 1994) demonstrate obesity and metabolic syndrome (Turek et al. 2005) as well as hypoinsulinemia and diabetes (Marcheva et al. 2010). These same dominant-negative mutant Clock-Δ19 mice also demonstrate an increase in cocaine reward, anxiety-related behavioral deficits, and behaviors similar to a human bipolar patient in a manic state (McClung et al. 2005; Roybal et al. 2007; Dzirasa et al. 2010, 2011). Multiple studies have suggested a link between CLOCK function and cancer (Alhopuro et al. 2010; Puram et al. 2016). In humans, a single nucleotide polymorphism (3111C) in the 3′ untranslated region (UTR) of CLOCK (Katzenberg et al. 1998) has been linked to difficulty losing weight (Garaulet et al. 2010), increased risk of insomnia (Serretti et al. 2003), and recurrence of major depressive episodes in bipolar disorder patients (Benedetti et al. 2003). In addition, a 4q12 copy number variation that includes CLOCK has been linked to autism spectrum disorder (ASD) (Sarachana et al. 2010; Griswold et al. 2012). Interestingly, variants in both the CLOCK direct transcriptional target PER1 and dimerization partner NPAS2 have also been associated with ASD (Nicholas et al. 2007; Neale et al. 2012). The preponderance of CLOCK involvement in these varied disease states suggests that this transcription factor plays a multitude of roles in regulating processes relevant to human health, some of which may be linked to circadian functions, and others that may be independent. Thus, understanding CLOCK function in regard to both circadian and noncircadian function requires further investigation for elucidating molecular mechanisms of many human disorders.

In order to investigate CLOCK function and the molecular machinery with which it interacts, most studies in mammals have focused on either nonneuronal tissues (such as the liver) or the core brain oscillator (the suprachiasmatic nucleus [SCN]) (Panda et al. 2002; Hughes et al. 2007, 2012; Koike et al. 2012). One recent study investigated CLOCK function in the neocortex and found that CLOCK plays an important role in regulating the critical period of visual cortex plasticity in mice (Kobayashi et al. 2015), supporting a critical role for CLOCK in the brain outside of the SCN. However, no study to date has specifically examined the role of CLOCK in the human neocortex. Such studies could lead to a deeper understanding of CLOCK-controlled targets associated with human brain disorders and could potentially open new therapeutic options for several neuropsychiatric diseases, such as sleep disorders, ASD, and epilepsy.

In line with this, previous work that identified CLOCK as a major hub gene in a human-specific cortical gene coexpression module also demonstrated significant enrichment for genes involved in neuropsychiatric disorders such as depression, schizophrenia, ASD, and seasonal affective disorder (Konopka et al. 2012a). Remarkably, the genes in this module do not significantly overlap with previously identified genes directly involved in circadian rhythms, suggesting identification of potentially novel targets of CLOCK specific to the human brain.

In line with this, it has been demonstrated that CLOCK expression is arrhythmic in the human neocortex (Li et al. 2013; Chen et al. 2016), supporting a potential noncircadian role for CLOCK in human cortical regions. Furthermore, the role of increased neocortical CLOCK expression as well as altered gene coexpression in human brain remains unknown. In this study, we determined the transcriptional networks regulated by CLOCK in the human brain through a combination of ChIP-seq and mRNA sequencing (RNA-seq). We identified several novel signaling cascades dependent on CLOCK expression, including those involved in metabolism, chromatin modifications, and neuronal migration. We functionally tested a number of these pathways and identified a novel role for CLOCK in regulating neuronal migration. Importantly, we demonstrated that CLOCK expression drives coexpression of genes that are involved in cognitive disorders and have human-specific patterns of expression. Together, these data provide new insights into human-specific patterns of gene expression and regulation in the brain, human brain evolution and specialization, and the transcriptional networks involved in human-specific cognitive disorders.

Results

Identification of CLOCK targets in the adult human neocortex

Direct transcriptional regulation of circadian genes by CLOCK has been described in multiple tissues in both mice and humans (Ripperger et al. 2000; Ueda et al. 2005; Ripperger and Schibler 2006; Alhopuro et al. 2010; Koike et al. 2012; Menet et al. 2012; Yoshitane et al. 2014; Puram et al. 2016). However, the transcriptional landscape governed by CLOCK in the human brain remains unexplored. To better understand the regulatory roles of CLOCK in this understudied environment, we profiled its DNA binding using ChIP-seq in the adult human neocortex (Supplemental Table S1).

Two neocortical regions in which CLOCK plays an important role in regulating the transcriptional landscape were chosen: Brodmann area 10 (BA10) and BA40 (Fig. 1A). BA10 (frontopolar prefrontal cortex) has special significance in regard to brain evolution, as it has been enlarged and structurally modified in humans throughout evolution (Semendeferi et al. 2011), is involved in higher-order cognitive tasks (Wendelken et al. 2011), and appears to be affected in neurological and psychiatric diseases, including ASD and schizophrenia (Dumontheil et al. 2008). Furthermore, BA10 is part of the frontal pole, and CLOCK was discovered to be significantly up-regulated and a major hub gene in a human-specific frontal pole module compared with other primates using weighted gene coexpression network analysis (WGCNA) (Konopka et al. 2012a). Moreover, network analysis of the human BrainSpan data set (http://developinghumanbrain.org) using weighted topological overlap (wTO) (Nowick et al. 2009) demonstrated that CLOCK has the most connections to other transcription factors in BA40 (inferior parietal cortex) relative to other cortical areas (Supplemental Fig. S1).

Figure 1.

Targets bound by endogenous CLOCK in the adult human brain. (A) CLOCK ChIP-seq was performed in two neocortical areas of the adult human brain: BA10 and BA40. (B) CLOCK ChIP-seq reads cluster near the transcription start sites (TSSs) of target genes and are highly enriched in the canonical CLOCK target E-box motif. (C) Most CLOCK ChIP-seq reads fall within ±1 kb of the TSS of target genes, although ∼10% are also found in distal intergenic regions. (D) Transcriptional targets of CLOCK in BA10 and BA40 significantly overlap (hypergeometric test, false discovery rate [FDR] = 2.7 × 10−85) but also include genes unique to each brain region. (E) Significant overlap of CLOCK ChIP-seq transcriptional targets from adult human brains with CLOCK ChIP-seq data from Puram et al. (2016). Hypergeometric test, FDR = 1.22 × 10−9. (F) Representative track files of BA10 and BA40 both show enriched CLOCK binding upstream of known canonical circadian targets (e.g., NR1D1) as well as noncircadian targets (e.g., ANKRD12). (G) Gene ontology analyses reveal that bound targets of CLOCK in either BA10 or BA40 are enriched for genes involved in RNA metabolism in addition to enrichment of genes important for canonical circadian functions.

ChIP antibody specificity was validated by performing ChIP-PCR in 293T cells after CRISPR-mediated CLOCK knockout (Mali et al. 2013; Wang et al. 2013). ChIP-PCR in wild-type 293T cells showed ∼16-fold enrichment of a known CLOCK-binding region over an adjacent site not transcriptionally occupied by CLOCK (Supplemental Fig. S2). This CLOCK ChIP-PCR enrichment was completely abolished in two different lines of CLOCK knockout 293T cells (Supplemental Fig. S2).

As expected for a transcription factor acting as a direct transcriptional activator, CLOCK ChIP reads clustered near the transcription start sites (TSSs) of target genes (Fig. 1B). Binding sites were enriched for the canonical E-box motif (P = 2.23 × 10−24 for BA10 and P = 1.54 × 10−29 for BA40, Fisher's exact test) CACGTG (Fig. 1B), as described previously for the CLOCK transcriptional complex (Ripperger and Schibler 2006).

Although most CLOCK ChIP reads fell within ±3 kb of the TSSs of target genes, a remarkable percentage of CLOCK ChIP reads were in distal intergenic regions (9.0% in BA10 and 11.4% in BA40) (Fig. 1B,C). These results suggest a potential regulatory role for CLOCK, extending beyond traditional binding to the promoter region immediately surrounding target genes. In addition, these results do not demonstrate a qualitative difference in CLOCK binding in BA10 compared with BA40 in terms of motif analysis or binding distribution.

We detected CLOCK binding upstream of 190 genes in BA10 and 159 genes in BA40, representing a highly significant overlap between cortical regions (65; hypergeometric test, false discovery rate [FDR] = 2.7 × 10−85) (Fig. 1D; Supplemental Table S1). These observations, along with the identification of the genes unique to each brain region, suggest that CLOCK could potentially play partially distinct but mostly overlapping roles in different neocortical areas in humans. BA40 has not been extensively compared from a gene expression perspective between humans and nonhuman primates, although the human cortical transcriptome is relatively homogenous outside of primary sensory areas (Hawrylycz et al. 2015). The similarity of CLOCK binding within these two cortical areas is therefore not particularly surprising, and any differences would need to be followed up in future studies.

To further validate our analysis of CLOCK-bound regions, we compared these CLOCK ChIP-seq data from adult human brains with previously generated lists of CLOCK targets from different tissues (Koike et al. 2012; Perelis et al. 2015; Puram et al. 2016). Our results demonstrated fewer called transcriptional targets (284 combined BA10 and BA40 CLOCK target genes) compared with two mouse Clock ChIP-seq studies from Perelis et al. (2015) in the mouse pancreatic β-cell line β-TC6 (1257 Clock target genes) and from Koike et al. (2012) in the mouse liver (3201 Clock target genes). Although we found fewer CLOCK target genes in human brains, the overlap with these two mouse Clock ChIP data sets was significant (85 genes in common with Koike et al. [2012] [FDR = 3.2 × 10−4] and 63 genes in common with Perelis et al. [2015] [FDR = 1.94 × 10−13]) (Supplemental Fig. S3). Puram et al. (2016) performed CLOCK ChIP-seq in human THP-1 cells (a monocytic cell line derived from an acute monocytic leukemia) and found the number of target genes (229) to be similar to that of the human brain data set, although their list of CLOCK target genes was based on the intersection of both CLOCK and ARNTL ChIP peaks (Puram et al. 2016). The 21 overlapping genes of these two data sets was also significant (hypergeometric test, FDR = 1.22 × 10−9) (Fig. 1E), providing additional confidence in the list of CLOCK targets identified in the human brain. Seven out of the 12 targets that overlap among all the ChIP-seq studies were canonical circadian targets (BHLHE41, CRY1, CRY2, DBP, NR1D1, PER1, and PER2), validating the ChIP-seq data sets and underscoring the idea that, while there is an essential role for CLOCK in regulating these few particular genes ubiquitously, the majority of CLOCK targets is tissue-specific.

Representative track files of BA10 and BA40 both showed strong CLOCK binding upstream of the known canonical targets NR1D1, HLF, and TEF (Fig. 1F; Supplemental Fig. S4A). Noncircadian CLOCK transcriptional targets such as ANKRD12, DNAJC16, and MAP2K7 were also strongly bound (Fig. 1F; Supplemental Fig. S4B). Intriguingly, NR1D1 has been associated recently with ASD and neurodevelopment (Goto et al. 2017). The PAR bZIP transcription factors HLF (Supplemental Fig. S4A), DBP (not shown), and TEF (Supplemental Fig. S4A) were strongly bound by CLOCK in the human brain, and dysregulation of both has been linked to epilepsy in mouse models (Gachon et al. 2004). A robust CLOCK-binding peak was shown in the promoter region of ANKRD12 (Fig. 1F), which was shown previously to coexpress with CLOCK specifically in the human frontal pole (Konopka et al. 2012a) and linked recently to childhood apraxia of speech (Peter et al. 2016). ChIP for overexpressed Flag-tagged CLOCK was used as a positive control in differentiated in vitro human neurons, with resulting data corroborating our in vivo human brain results (Fig. 1F).

Gene ontology (GO) analysis was performed to identify functional categories of all of the combined CLOCK ChIP targets from both BA10 and BA40 (Supplemental Table S1). Due to high similarity (sim = 0.378) (Schlicker et al. 2006) between the GO lists resulting from the ChIP targets of BA10 and BA40, we integrated both gene lists into a consensus CLOCK ChIP target gene set (ConsensusChIPTargets) for downstream analyses (Supplemental Fig. S5). Interestingly, these neocortical targets of CLOCK were enriched in genes involved in RNA metabolism in addition to an enrichment of genes important for canonical circadian function (Fig. 1G). These experiments are the first to identify bound CLOCK targets in the human neocortex and leave the door open for further investigation of the functional effects of CLOCK control of these direct transcriptional targets.

Identification of DEGs following CLOCK knockdown in human neurons

To test the underlying mechanism of CLOCK regulation of gene expression in a genetically manipulable human neuronal environment, we performed RNA-seq in differentiated human neuronal progenitors (hNPs) after CLOCK knockdown and analyzed global gene expression changes (Supplemental Fig. S6A–C; Supplemental Table S2).

CLOCK is a canonical circadian transcription factor shown to have rhythmic transcriptional activity in some tissues (Koike et al. 2012), although human gene expression studies have suggested that CLOCK itself does not oscillate on the RNA level in the human neocortex (Li et al. 2013; Chen et al. 2016). However, due to the scarcity of circadian gene expression studies in the human brain and the abundance of studies showing rhythmicity of CLOCK mRNA transcripts in different tissues (Gekakis et al. 1998; Steeves et al. 1999; Panda et al. 2002; Ueda et al. 2005; Rath et al. 2012), we performed a full 2-d time course, taking samples every 4 h after synchronization of the neurons with dexamethasone (McNamara et al. 2001) to capture any potential circadian phase-dependent gene expression changes (Fig. 2A). We compared CLOCK knockdown cells with cells expressing a scrambled nontargeting shRNA. In both cases, cells were treated equivalently with dexamethasone and compared only with each other, negating the need for an “untreated” group.

Figure 2.

DEGs in differentiated human neurons following CLOCK knockdown. (A) RNA-seq was performed on differentiated human neuronal samples collected every 4 h for 2 d following CLOCK knockdown. (B) A heat map of DEGs demonstrates robust clustering by genotype rather than time point. (C) An MA plot showing all DEGs with log2 (fold change) > |0.3|, with up-regulated DEGs in blue (379 in total) and down-regulated DEGs in red (185 in total). All DEGs meeting this criterion have an FDR < 0.05. The bottom left inset shows the distribution of genes by FDR, with the dashed blue line representing the FDR cutoff of <0.05. The bottom right inset box plot displays the average CLOCK expression in control and CLOCK knockdown samples, representing a log2 fold change of −1.502 in knockdown samples relative to controls. (D) Down-regulated DEGs significantly overlap with CLOCK ChIP-seq targets in the human brain (hypergeometric test, FDR = 2.07 × 10−3) and in human neurons (hypergeometric test, FDR = 1.40 × 10−6). (E) Up-regulated DEGs are significantly enriched for genes involved in cell migration and carbohydrate metabolism.

We found 564 DEGs (185 down-regulated and 379 up-regulated) after CLOCK knockdown (Fig. 2B,C). DEGs were defined using an FDR ≤ 0.05 and log2 (fold change) ≥ |0.3| as done previously (Araujo et al. 2015), and permutation testing supported these filters (see the Materials and Methods; Supplemental Material). Interestingly, these DEGs displayed robust clustering by genotype rather than time point, and neither cluster of DEGs oscillated in expression relative to a 24-h circadian cycle (Fig. 2B). Furthermore, common algorithms used to call cycling genes, such as ARSER and JTK_CYCLE (Hughes et al. 2010; Yang and Su 2010; Wu et al. 2016), did not detect any cycling genes in the entire data set (including CLOCK and ARNTL in the control cells), in agreement with quantitative RT–PCR (qRT–PCR) for a small number of genes (data not shown), suggesting that the DEGs regulated by CLOCK display noncycling expression specifically in human neurons as hypothesized. However, the differentiated neurons did display some rhythmic expression of an ARNTL promoter driving luciferase in both control and CLOCK knockdown genotypes, albeit at extremely low amplitude (Supplemental Fig. S7A,B). These differences likely reflect the continuous repeated lumicycle measurements in the same samples versus regular, but relatively infrequent, RNA-seq sampling and the difference in single-gene sensitivity in luciferase screening compared with RNA-seq. To further characterize the expression patterns of some of the core circadian factors, we also examined differential expression of these factors at each time point in our data set (Supplemental Fig. S8). We found that with CLOCK knockdown, NR1D1, NR1D2, PER1, HLF, and TEF tend to be down-regulated, supporting positive regulation of these genes by CLOCK. Beyond that conclusion, we cannot speculate further, as many of these genes regulate one another through positive and negative feedback loops. However, it is clear that reduction of CLOCK itself is sufficient to disrupt expression of several core factors.

All differential expression and coexpression network analyses were derived from RNA-seq data from a single hNP cell line differentiated into neurons and using one specific shRNA to knock down CLOCK. We further validated our RNA-seq results by reproducing similar expression changes in representative genes using both an alternative shRNA to CLOCK and an independent hNP line (Supplemental Fig. S9). We found that many of the same gene expression changes that occur following CLOCK knockdown in this different hNP cell line differentiated into neurons as well as with a distinct shRNA that targets a separate region of CLOCK (i.e., up-regulated DEGs CTGF and PRDX3). Differential expression of some canonical targets, such as NR1D1, was recapitulated upon knockdown of CLOCK heterodimer partner ARNTL (Supplemental Fig. S9). However, these differentially expressed results were not preserved by knocking down either ARNTL or the CLOCK paralog NPAS2 (i.e., PRDX3) (Supplemental Fig. S9), providing further evidence of the potential noncircadian role of CLOCK in contrast to the other core circadian factors and the overall specificity of the results.

Because CLOCK is generally considered to be an activator of transcriptional targets, we predicted significant overlap between down-regulated DEGs following CLOCK knockdown and the in vivo ConsensusChIPTargets from the brain as well as the in vitro ChIP-seq targets from the same neuronal cell line. As expected, the overlap was significant between down-regulated DEGs and both the in vivo brain and in vitro neuron CLOCK ChIP targets (hypergeometric test, FDR = 2.07 × 10−3 and FDR = 1.40 × 10−6, respectively) (Fig. 2D). There were 10 genes in common between human brain CLOCK ChIP targets and down-regulated DEGs (BHLHE41, DBP, DNAJC16, HLF, MAP2K7, NR1D1, NR1D2, PER1, PER3, TEF), providing a high-confidence list of CLOCK transcriptional target genes with both in vivo and in vitro relevance (Figs. 1F, 2D; Supplemental Fig. S4A,B). ChIP for overexpressed Flag-tagged CLOCK significantly overlapped with the gene target list from both the human brain (hypergeometric test, FDR = 1.13 × 10−23) and the down-regulated DEGs (1.40 × 10−6) (Fig. 2D). DBP, HLF, and TEF, as mentioned previously, have been linked to regulation of neuronal excitability in the brain (Gachon et al. 2004). Interestingly, MAP2K7 expression was decreased in the prefrontal cortex of human patients with schizophrenia, and mice with only one copy of Map2k7 exhibited impaired working memory reminiscent of a schizophrenia-like cognitive phenotype (Winchester et al. 2012; Openshaw et al. 2017). Finally, DNAJC16, a poorly studied member of the heat-shock protein family, was also a transcriptional target of CLOCK in the human brain (Supplemental Fig. S4B) and was shown to be dysregulated in human neurons following CLOCK misexpression (Fig. 2D).

In contrast, up-regulated DEGs showed no such enrichment in ConsensusChIPTargets (hypergeometric test, FDR = 0.9752), suggesting that CLOCK is influencing expression of these genes indirectly through intermediate transcriptional networks. The strong down-regulation of NR1D1 (itself a potent transcriptional repressor) (Harding and Lazar 1995; Preitner et al. 2002; Ko and Takahashi 2006) following CLOCK knockdown could provide a partial explanation for the greater number of up-regulated DEGs compared with down-regulated DEGs. Functionally, up-regulated DEGs were enriched in cell localization and cell projection function, representing potentially novel groups of genes impacted by CLOCK expression (Fig. 2E; Supplemental Table S2). Up-regulated DEGs were also functionally enriched in metabolic processes (Fig. 2E), potentially echoing the functional role of ConsensusChIPTargets in regulating RNA metabolism (Fig. 1G; Supplemental Table S2). Altogether, these observations suggest that CLOCK might have specific noncycling expression in human neurons, with a crucial noncircadian role in regulating genes implicated in cell migration and neuronal projections.

CLOCK expression is related to coexpression of genes involved in epigenetic modifications

To determine the network of genes regulated by CLOCK as well as prioritize the identified DEGs, we applied WGCNA (Langfelder and Horvath 2008) to the RNA-seq data set. We identified 21 modules containing highly coexpressed genes in total (Supplemental Table S3). Three of these modules were related to CLOCK expression or enriched for genes involved in neuropsychiatric diseases such as ASD and intellectual disability (Fig. 3).

Figure 3.

WGCNA reveals modules of coexpressed genes related to CLOCK expression in human neurons. Heat map plotting significant enrichment of gene sets related to CLOCK expression or human neuropsychiatric disorders in CLOCK modules (CMs) 1–3.

CLOCK module 1 (CM1) contained CLOCK itself and was highly enriched in down-regulated DEGs, and, as such, the module eigengene associated negatively with CLOCK knockdown across all time points tested (Figs. 3, 4A,B). Although CM1 was not significantly enriched for the neocortical ConsensusChIPTargets (hypergeometric test, FDR = 0.1284), some of the most highly connected genes (e.g., DNAJC16, BTBD7, and NR1D2) were among the ConsensusChIPTargets (Fig. 4B). CM1 was enriched in genes involved in chromosome organization (Supplemental Table S3), which, in combination with the reported histone acetyltransferase (HAT) function of CLOCK itself (Doi et al. 2006), suggested that CLOCK may play a role in regulating epigenetic modifications (Fig. 4C). To test this hypothesis, potential changes in the open chromatin state following CLOCK knockdown in differentiated human neurons were assayed genome-wide by assay for transposase-accessible chromatin (ATAC) sequencing (ATAC-seq). Surprisingly, we found few significant differentially open regions (Fig. 4D), although positive controls for the assay, such as examining the CLOCK gene with overexpression of CLOCK, demonstrate its functionality (data not shown). The lack of alteration in the open chromatin state with CLOCK knockdown applied to both the genome-wide data and the regulatory regions associated with genes in CM1. These data suggest that CLOCK expression does not directly affect the open chromatin state on the time frame assayed, and the biological implications of the coexpressed genes in CM1 remain to be determined.

Figure 4.

Identification of CLOCK-related gene expression involved in epigenetic modifications. (A) The module eigengene of CM1 correlates negatively with CLOCK knockdown across all time points tested. (B) CM1 contains CLOCK itself and is significantly enriched in down-regulated DEGs. (C) CM1 is significantly enriched for genes involved in chromosome organization and ubiquitin pathway modification. (D) MA plot of ATAC-seq binding affinity comparing control and CLOCK knockdown samples that detected few differentially open chromatin regions with an FDR < 0.05.

CLOCK knockdown results in up-regulation of genes involved in neuronal migration and increases neuronal migratory distance in vitro

In contrast to CM1, CM2 was highly enriched in up-regulated DEGs, and the aggregate expression of the genes in this module associated positively with CLOCK knockdown (Figs. 3, 5A,B). The genes in CM2 were involved in cell migration and neurogenesis (Fig. 5C; Supplemental Table S3), similar to the functional analysis results from up-regulated DEGs overall (Fig. 2E; Supplemental Table S2).

Figure 5.

CLOCK expression regulates neuronal migration. (A) The module eigengene of CM2 correlates positively with CLOCK knockdown. (B) CM2 is significantly enriched in up-regulated DEGs. (C) CM2 genes are significantly enriched for genes involved in cell migration and neurogenesis. (D) Human neural progenitors transduced with control virus (RFP+) form spherical “neurospheres” and migrate away from the core on laminin-coated plates. All cells that migrated past the white dashed line, representing 200% of the inner core size, were counted in this analysis. n = 350–500 cells per neurosphere. Red dashed lines indicate measurement from the inner core to the outermost control virus-positive cell (for clarity, not all measurements are shown). Bar, 200 µm. (E) Neurospheres comprised of cells transduced with CLOCK knockdown virus migrated further from the core than did neurospheres made up of cells infected with control virus. Two-tailed Mann-Whitney test, P = 3.79 × 10−2. (F) Single neurospheres were created by mixing equal numbers of cells infected separately with control virus (GFP+) and CLOCK knockdown virus (RFP+). All cells that migrated past the white dashed line, representing 200% of the inner core size, were counted in this analysis. n = 350–500 cells per neurosphere. Green dashed lines show the distance between the inner core and the outermost control virus-infected cell, and red dashed lines represent measurement from the core to the outermost CLOCK knockdown virus transduced cell (for clarity, not all measurements are shown). Bar, 200 µm. (G) The same phenotype of greater migratory distance traveled with CLOCK knockdown relative to control was observed with greater significance using “within neurosphere replicate” controls. Two-tailed paired t-test, P = 3.0 × 10−3.

Several specific genes in CM2 that were up-regulated following CLOCK knockdown could be playing key roles in regulating neuronal migration. PRDX3 is the second most interconnected gene in CM2, was confirmed as up-regulated with an independent hNP line and alternative shRNA to CLOCK (Supplemental Fig. S9), and has been shown to mediate cell migration in the HepG2 hepatocellular carcinoma cell line (Liu et al. 2016). DISC1 is heavily implicated in both the schizophrenia and ASD disease process, and a mutant form of this gene has been linked to abnormal neocortical development due to migration defects (Kamiya et al. 2005). MET is also associated with both ASD and defects in neuronal migration in both animal models and post-mortem human studies (Campbell et al. 2007; Elsen et al. 2009; Mukamel et al. 2011). Finally, a comparative genomics study examining mouse and human gene expression in radial glia (RGs) found higher expression of the CM2 gene PDGFRβ (and its ligand, PDGFD) in the human brain, and pharmacological inhibition of PDGFD–PDGFRβ signaling prevented normal cell migration of neocortical RGs in human, but not mouse, slice cultures (Lui et al. 2014). Together, these data suggested that CLOCK regulates the expression of genes related to migratory function in these human neuronal cells.

To test the functional significance of the cell migration/localization GO association of DEGs and CM2 genes (Figs. 2E, 5C), we performed a neuronal migration assay. Neurospheres of human neural progenitor cells transduced with either control or CLOCK knockdown lentivirus formed into a spherical single unit of cells, and cells were permitted to migrate away from the core over 2 d (Fig. 5D). Neurospheres comprised of cells transduced with CLOCK knockdown virus migrated further from the core than did neurospheres made up of cells infected with control virus (Fig. 5E). However, each individual neurosphere varied significantly in core size and average migration distance, potentially adding additional variance beyond the manipulation of the desired experimental variable. Single neurospheres were therefore created by mixing equal numbers of cells transduced separately with either control virus (GFP+) or CLOCK knockdown virus (RFP+) (Fig. 5F). This approach generated within-replicate controls and revealed a more robust phenotype statistically (Fig. 5G). Furthermore, this additional evidence suggested that a secreted factor regulating migratory distance is a less likely explanation for these results. We further performed an additional independent assay of neuronal migration following CLOCK knockdown using the established transwell assay (Kaplan et al. 2005; Ma et al. 2007). These results corroborated our neurosphere data and further supported a role of CLOCK in regulating neuronal migration (Supplemental Fig. S10).

CLOCK-relevant human-specific hub genes are enriched in genes implicated in neurodevelopmental disorders

To identify the relevance of the identified CLOCK target genes to human brain evolution, we first examined whether there was significant overlap of CLOCK target genes and previously identified lists of genes with human-specific expression (Caceres et al. 2003; Khaitovich et al. 2004, 2005; Babbitt et al. 2010; Konopka et al. 2012a; Liu et al. 2012). However, we found no significant enrichment of genes with human-specific expression in ConsensusChIPTargets or DEGs following CLOCK knockdown (hypergeometric test, FDR = 0.9989 and FDR = 0.7904, respectively).

Because cognitive disorders are thought to have arisen partly as a consequence of human brain evolution (Fontenot and Konopka 2014), we next examined whether any CMs were enriched in genes involved in cognitive disorders. CM3, a module with aggregate expression of the genes associating positively with CLOCK knockdown (Fig. 6A), was enriched in high-confidence “scored” ASD-associated genes as well as genes associated with intellectual disability (Figs. 3, 6B). CM3 was enriched in genes involved in neuronal development and projection assembly in addition to genes implicated in several neuropsychiatric disorders as mentioned above, such as ASD, intellectual disability, and epilepsy (Fig. 6C; Supplemental Table S3). In particular, the CM3 hub gene SRGAP3 has been associated specifically with X-linked intellectual disability and is involved in the Slit–Robo pathway regulating axonal branching and neuronal migration (Endris et al. 2002).

Figure 6.

CM3 is enriched in genes relevant to human neuropsychiatric health. (A) The expression of CM3 genes is positively associated with CLOCK knockdown. (B) CM3 is significantly enriched in up-regulated DEGs and disease-related genes. (C) CM3 is significantly enriched for genes involved in a number of neuropsychiatric disorders, such as ASD, intellectual disability, and epilepsy. (D) The most interconnected genes (or hubs) of CM3 are significantly enriched for genes with human-specific patterns of brain expression compared with other primates. “CM3 genes,” bidirectional Wilcoxon signed-rank test, P = 1.29 × 10−3. This result is specific to CM3, as the hub genes in the remaining modules do not demonstrate this enrichment. “All WGCNA genes,” bidirectional Wilcoxon signed-rank test, P = 1.99 × 10−1.

None of the modules from the WGCNAs was significantly enriched for either human-specific genes or the genes from the original human frontal cortex module that identified CLOCK as having human-specific connectivity (Fig. 3; Konopka et al. 2012a). However, the techniques and sample sizes are quite varied for the identification of human-specific genes, likely making these lists of genes underpowered and/or noisy (see the Materials and Methods; Supplemental Material). We therefore focused on the most interconnected genes of the modules (“hub” genes), as these genes drive the coexpression of the other genes within a given module (Fontenot and Konopka 2014). We found that the hub genes of CM3 tended to have human-specific patterns of expression compared with other primates, which was a feature unique to this module among all of the modules in our data set (Fig. 6B,D). As CLOCK itself showed a human-specific pattern of expression in the human cortex (Babbitt et al. 2010; Konopka et al. 2012a), these results support the hypothesis that CLOCK plays a key role in regulating gene expression networks relevant to human cognition and related disorders.

Discussion

To identify human-specific neuronal targets, one is limited by the availability of cells with in vivo relevance. We used a well-characterized system of hNPs that has been shown previously to recapitulate in vivo brain development (Konopka et al. 2012b; Stein et al. 2014). These cells are most closely related to forebrain neurons and therefore are more relevant to cortical neurons than to brain regions with oscillatory functions such as the SCN. Using this human cellular system, we can directly manipulate CLOCK expression—something that cannot be done in vivo with humans. We additionally used post-mortem adult human brain tissue in addition to the same in vitro neuronal model to carry out CLOCK ChIP-seq to address the in vivo relevance of our results. The overlap of these ChIP-seq data together with the neuronal RNA-seq data provides some validation of directly bound transcriptionally regulated targets, while the ChIP-seq targets on their own provide an additional set of genes for further examination in cellular or other model systems. The overlap of these data at both the individual gene level and the functional pathway level is even more remarkable, as we were comparing data from cells derived from developing human brains with data from adult human post-mortem tissue. Limitations regarding the availability of sufficient quantities of human fetal brains for ChIP-seq or of primary adult human neurons that can be cultured make carrying out comparable time points for in vivo and in vitro studies formidable. Using available human fetal brain gene expression data sets (Miller et al. 2014), we found that CLOCK is highly expressed throughout human fetal brain development (data not shown), suggesting that it could be playing a role in human cortical development.

Determining the human or brain specificity of these results is also challenging, as there have been few studies that have examined genome-wide direct targets of CLOCK with which to compare. Although there were significant overlaps, the CLOCK ChIP-seq data from adult human brains encompass a mostly unique list of transcriptional targets compared with the CLOCK ChIP-seq studies carried out in mouse livers (Koike et al. 2012), mouse pancreatic β cells (Perelis et al. 2015), and human THP-1 cells derived from an acute monocytic leukemia patient (Puram et al. 2016). We believe that most of these differences are largely due to the choice of organs studied, species, and antibodies. However, there are also important technical differences in the way the ChIP experiments were carried out. Most notably, our data are derived from samples single-cross-linked with formaldehyde, while Koike et al. (2012), Perelis et al. (2015), and Puram et al. (2016) all double-cross-linked samples with either disuccinimidyl glutarate (DSS) and formaldehyde or ethylene glycol bis (succinimidyl succinate [EGS]) and formaldehyde. DSS and EGS can cross-link protein–protein interactions, and thus their use in ChIP-seq protocols could identify both direct DNA binding and indirect binding through the action of CLOCK as a cofactor (Zeng et al. 2006). Furthermore, as the genes expressed vary across these tissues and species, it is not surprising that regulation of gene expression would also vary.

Both CLOCK ChIP results in human brains and differentially expressed analysis following CLOCK knockdown in human neurons suggest a potential functional role for CLOCK in regulating metabolic processes in RNA (Fig. 1G) and carbohydrate (Fig. 2E) metabolism, respectively. This regulatory role is especially interesting from an evolutionary perspective, as many studies have suggested that change in metabolic regulation was necessary for the development of the enlarged human cortex relative to other primates (Aiello and Wheeler 1995; Uddin et al. 2004; Bauernfeind et al. 2014; Kuzawa et al. 2014). Future studies could measure the metabolic products or rate following CLOCK manipulation in in vitro human systems or mouse models misexpressing CLOCK. Studies in nonhuman primate cells or tissues would also be informative, but, unfortunately, access to the comparable materials is extremely limited.

Cyclical transcriptional fluctuations have been shown to correlate with epigenomic changes in chromatin structure over circadian time (Masri and Sassone-Corsi 2010). Previous reports have suggested that CLOCK demonstrates intrinsic HAT activity and can acetylate histone H3 (Doi et al. 2006; Hirayama et al. 2007). Furthermore, CLOCK and ARNTL interact with coactivators p300 and CREB-binding protein (CBP) to acetylate histones and prepare an accessible chromatin state for transcription (Etchegaray et al. 2003; Curtis et al. 2004; Lee et al. 2008; Hosoda et al. 2009). We therefore hypothesized that any effects of increased CLOCK expression in the human lineage could perhaps be explained in part by HAT function. These ideas were extended by the GO enrichment of CM1 genes in categories relevant to epigenetic modification (Fig. 4C). We were therefore surprised by the negative results obtained when we tested this hypothesis using ATAC-seq after CLOCK knockdown (Fig. 4D). However, these negative results leave open the possibility that some yet to be determined epigenetic modifications could be occurring despite the stability of open chromatin regions following CLOCK knockdown. These untested modifications potentially performed by CLOCK target genes include direct DNA methylation (Harris et al. 2010) as well as phosphorylation, ubiquitination, and sumoylation of histone proteins (Geiss-Friedlander and Melchior 2007; Bannister and Kouzarides 2011). Modifications such as these could potentially have been missed by our ATAC-seq experiments and warrant further study.

Mature neurons in the human brain must migrate far from their origins to reach their programmed positions in the adult neocortex, and defects in this process can result in many neurological diseases (Gleeson and Walsh 2000). Interestingly, CLOCK is a member of the basic helix–loop–helix (bHLH) family of transcription factors, and this family has been linked to regulation of cell migration in previous studies (Ge et al. 2006). Furthermore, another core circadian gene and direct CLOCK transcriptional target (Fig. 1F), NR1D1, was shown to regulate neuronal migration in the cerebral cortices of mice, and mutations in this gene were linked to ASD in human patients (Goto et al. 2017). Upstream CLOCK control of NR1D1 could potentially play a role in this corroborating migratory phenotype. Such a result would be in line with up-regulated genes after CLOCK knockdown being associated with neuronal migration via dysregulation of NR1D1. Indeed, repeating our neurosphere assay of migration with knockdown of NR1D1 showed that, as predicted, NR1D1 knockdown promoted further neuronal migration in a fashion similar to that of CLOCK knockdown experiments (Supplemental Fig. S11). The genes involved in regulation of neuronal migration may not all be direct targets of CLOCK, but understanding the networks of gene expression downstream from CLOCK in human neurons in particular significantly contributes to our knowledge of the molecular mechanism of human brain development and evolution.

Our study is the first to report regulation of neuronal migration by CLOCK. Interestingly, a recent study suggested that CLOCK plays a vital role in regulating the critical period of mouse visual cortex plasticity (Kobayashi et al. 2015) but did not posit a specific mechanism. A change in neuronal migration regulation downstream from CLOCK could potentially provide such an explanation by altering the timing and position of the normal developmental landscape and could be investigated further in mouse models with altered CLOCK levels. In contrast, it is possible that the negative regulation of neuronal migration by CLOCK could be relevant only to human brain development. There are significant differences in the human brain, in particular the increased surface area of the human neocortex. Thus, the tradeoff of migration versus neural stem cell division plays an important role in directing human cortical expansion (Geschwind and Rakic 2013). Based on the expression of CLOCK in the developing human brain, it is tempting to speculate that increased CLOCK expression is important for regulating molecular pathways involved in maintaining a proliferative state rather than promoting migration. Thus, while our data are provocative and the first of their kind to link CLOCK to human brain development, many compelling avenues of inquiry remain for future studies.

Elucidating human-relevant transcriptional cascades in the brain should provide insight into both human brain evolution and disease. To our knowledge, the study detailed here is the first of its kind to identify the transcriptional program downstream from CLOCK in human neurons. These findings provide context to the previous observation of CLOCK demonstrating increased expression in the human neocortex, supporting a noncircadian role for CLOCK in this region. We found that CLOCK regulates the expression of several functional categories of genes outside of the known circadian circuits, most notably genes involved in metabolism and cell migration. We cannot rule out that although CLOCK itself is not robustly cycling, it could be functioning in a circadian manner through cyclical interaction with its dimerization partner, ARNTL, which is the top-ranked cycling gene in several human brain regions, including neocortex (Li et al. 2013). Such a hypothesis is further supported by our identification of the enrichment of an E-box motif in the ChIP-seq data as well as the coexpression of CLOCK with two other proteins that robustly cycle in the human brain: PER2 and NR1D2 (Supplemental Fig. S1). Future studies that assess CLOCK protein interactions in the human brain throughout circadian time should address these possibilities. We provided additional experimental evidence to demonstrate a repressive role for CLOCK in regulating human neuronal migration. Finally, we also uncovered an enrichment of genes associated with cognitive disorders or human-specific brain expression patterns that are related to CLOCK expression in human neurons. Together, these data support an evolutionarily relevant role for CLOCK to direct signaling cascades important for brain development that are at risk in cognitive diseases.

Materials and methods

CLOCK ChIP

ChIP experiments were carried out similar to those in Koike et al. (2012). Please see the Supplemental Material for more information.

hNPs and human neurons

hNPs were purchased from Lonza and cultured as described previously (Konopka et al. 2012b). Please see the Supplemental Material for more information.

ChIP-seq and RNA-seq library construction and sequencing

ChIP-seq and RNA-seq libraries were prepared in-house as described previously (Takahashi et al. 2015). For RNA-seq libraries, mRNA was isolated from randomized samples using poly-A selection, and strand-specific libraries were generated. Sequencing was performed on randomly pooled samples by the McDermott Sequencing Core at University of Texas Southwestern on an Illumina NextSeq 500 sequencer. Single-end 75-base-pair (bp) reads were generated.

Differential expression analyses

Differential expression between control and CLOCK knockdown samples was assessed using the DESeq2 package in R (Supplemental Table S2; Anders and Huber 2010). The expression matrix contained the 12,183 protein-coding genes that passed the RPKM (reads per kilobase per million mapped reads) cutoff as described in the Supplemental Material. Linear regression was performed using DESeq2 to remove covariate variables: RNA integrity number (RIN) value and library construction batch. All DEGs with an FDR ≤ 0.05 and log2 (fold change) ≥ |0.3| (Araujo et al. 2015) were retained. A permutation test was also applied using 1000 permuted experiments. None of these permuted analyses showed the same genes differentially expressed (permutation P < 0.001).

Coexpression network analyses

To identify modules of coexpressed genes in the RNA-seq data, we carried out WGCNA (Supplemental Table S3; Langfelder and Horvath 2008). RPKM values were filtered as described in the Supplemental Material, and log2 (RPKM + 1) was used as input data. Please see the Supplemental Material for more information.

Statistical analysis and code availability

Statistical analyses were performed using R. All code is available on request.

Please see the Supplemental Material for details of all other methods.

All work using deidentified human materials have been granted an exemption from requiring ethics approval from University of Texas Southwestern Medical Center.

Availability of data and material

The NCBI Gene Expression Omnibus (GEO) accession number for the ChIP-seq, RNA-seq, and ATAC-seq data reported in this manuscript is GSE96659.