Potential contribution of tandem circadian enhancers to nonlinear oscillations in clock gene expression.

Limit-cycle oscillations require the presence of nonlinear processes. Although mathematical studies have long suggested that multiple nonlinear processes are required for autonomous circadian oscillation in clock gene expression, the underlying mechanism remains controversial. Here we show experimentally that cell-autonomous circadian transcription of a mammalian clock gene requires a functionally interdependent tandem E-box motif; the lack of either of the two E-boxes results in arrhythmic transcription. Although previous studies indicated the role of the tandem motifs in increasing circadian amplitude, enhancing amplitude does not explain the mechanism for limit-cycle oscillations in transcription. In this study, mathematical analysis suggests that the interdependent behavior of enhancer elements including not only E-boxes but also ROR response elements might contribute to limit-cycle oscillations by increasing transcriptional nonlinearity. As expected, introduction of the interdependence of circadian enhancer elements into mathematical models resulted in autonomous transcriptional oscillation with low Hill coefficients. Together these findings suggest that interdependent tandem enhancer motifs on multiple clock genes might cooperatively enhance nonlinearity in the whole circadian feedback system, which would lead to limit-cycle oscillations in clock gene expression.

INTRODUCTION

Circadian rhythms are essential biological processes in almost all organisms, ranging from unicellular to multicellular ones (Young and Kay, 2001). In mammals, chronic desynchronization between physiological and environmental rhythms carries a significant risk of diverse disorders, ranging from sleep disorders to diabetes, cardiovascular diseases, and cancer (Wijnen and Young, 2006; Bass, 2012; Sahar and Sassone-Corsi, 2012). The molecular machinery for circadian oscillation is cell autonomous (Rosbash, 1998) and ubiquitous throughout the body (Hastings ). Peripheral clocks enable autonomous circadian gene expression in a tissue-specific phase and manner. Just like radio-controlled clocks, the suprachiasmatic nucleus functions as a standard clock and then adjusts peripheral clocks to prevent internal desynchronization (Schibler and Sassone-Corsi, 2002). The circadian clockwork consists of a negative feedback loop of transcription, which drives cell-autonomous circadian oscillation in clock gene expression (Dunlap, 1999). In mammals, the BMAL1 and CLOCK complex activates transcription of the clock and clock-related genes via E-box elements. Subsequently Period (PER), together with Cryptochrome (CRY), functions to negatively regulate this complex (Reppert and Weaver, 2001; Ko and Takahashi, 2006). A limit cycle—a hallmark of self-sustained oscillations generated from nonlinear systems—provides a theoretical basis for cell-autonomous circadian gene expression (Goodwin, 1965). A sigmoidal characteristic in the transcription of clock genes has been considered a primordial nonlinearity for the realization of autonomous transcriptional oscillation. The underlying mechanism of this nonlinear transcription, however, remains elusive.

As shown in several biological processes (Ferrell, 1996), it is possible that transcriptional regulation of clock genes also obeys the Hill equation. Circadian transcription-translation feedback loops are therefore generally simulated in silico using the Goodwin model, which is based on the Hill function (Goodwin, 1965). However, realizing the degree of nonlinearity considered essential for autonomous transcriptional oscillation is presumed to require steep transcriptional activity kinetics, which in turn requires upward adjustment of the Hill coefficient to experimentally unrealistic values. Although this difficulty can be avoided by mathematically seeking an alternative source of nonlinearity in biological processes, such as multistep phosphorylation and nuclear translocation (Gonze and Abou-Jaoude, 2013), evidence for the presence of sufficient nonlinearity in these processes remains to be provided. Thus a deeper understanding of biological processes that contribute to the enhancement of nonlinearity requires further detailed analyses of the molecular mechanism of circadian transcription.

Although CRY2 is a core component of the mammalian circadian clock (van der Horst ), the mechanism of the circadian transcription of this gene has yet to be analyzed in detail. In the present study, we find that cell-autonomous circadian transcription of the Cry2 gene requires a functionally interdependent tandem E-box motif. Previous studies found tandem motifs of E-box commonly in the transcription-regulatory region of several other clock genes and indicated the role of these motifs in enhancing circadian amplitude (Akashi and Takumi, 2005; Nakahata ; Paquet ; Rey ). However, enhancing amplitude does not provide evidence for limit-cycle generation. In this study, we indicate that the interdependent behavior of circadian enhancer elements including not only E-boxes but also retinoic acid–related orphan receptor (ROR) response elements (ROREs) might contribute to the limit-cycle oscillation by increasing transcriptional nonlinearity.

RESULTS

To identify components necessary for cell-autonomous circadian transcription of the hCry2 gene, we isolated a 3-kbp 5’-flanking region upstream from the hCry2 transcription start site, subcloned it into the upstream of the luciferase gene, and monitored the hCry2-driven luciferase activity in real time using a cell-based system (Figure 1A). To exclude systemic factors such as blood-borne factors and body temperature, it was necessary to examine cell-autonomous transcriptional regulatory mechanisms under an in vitro culture condition. The data clearly showed a robust circadian oscillation of the bioluminescence in U2OS human osteosarcoma cells (Figure 1B), indicating that the cloned genomic region contained elements essential for circadian transcription. To narrow the transcriptional regulatory region required for the cell-autonomous circadian transcription of hCry2, we used a deletion series of hCry2 driven-reporter constructs and monitored transcriptional fluctuation in real time (Figure 1B). Transfected cells showed almost identical phases and amplitudes for all deletion constructs. The results indicated that the region spanning from −249 to +41 is necessary and sufficient for cell-autonomous circadian transcription of the hCry2 gene and that the region contains regulatory elements for it. To further narrow the genomic region essential for cell-autonomous circadian transcription, we constructed a deletion series of the hCry2 (–249 to +41) construct (Figure 1C). The results strongly suggested that the essential regulatory elements are contained in the region from –175 to –116 (60–base pair length; Figure 1D). The data from experiments with the deletion constructs indicated that the 60–base pair region is functionally divided into two—a 43–base pair region (from −158 to −116) contains elements that enhance oscillation amplitude, and a 17–base pair region (from −175 to −159) is indispensable for the generation of transcriptional oscillation. To identify consensus sequences for transcription factor binding within the latter region, we performed an in silico search by sequence alignment and found interspecies-conserved E-boxes within the 17–base pair region (E1 and E2; Figure 1C). Of interest, these two E-boxes are tandemly repeated and spaced by five nucleotides. Although a previous comprehensive bioinformatics analysis of evolutionarily conserved genomic regions demonstrated the presence of clock-controlled elements on clock genes, no elements, such as E-boxes, that might potentially control Cry2 transcriptional oscillation were identified (Ueda ).

FIGURE 1:

Mapping of the upstream region of hCry2 required for transcriptional oscillation. (A) Schematic representation of deletion mutants of a 3-kbp 5’-flanking region upstream from the hCry2 transcription start site. Here +1 corresponds to the transcription start site. (B) Cell culture–based luminescent monitoring was performed with the hCry2 constructs indicated. U2OS human osteosarcoma cells were transfected with hCry2-luc constructs and then stimulated with dexamethasone. After the synchronization, in the presence of luciferin, light emission was measured and integrated for 1 min at intervals of 15 min. The maximum bioluminescence was set to 0.6. The warm and cold colored lines represent data from three independent biological replicates for the full-length hCry2-luc (−2925 to +41) construct and its deletion mutants, respectively. Data sets obtained from three experiments were detrended by subtracting the 24-h running average from raw data. (C) Schematic representation of a deletion series of the hCry2-luc (−249 to +41) construct. E1 and E2 represent putative E-boxes. (D) U2OS human osteosarcoma cells were transfected with hCry2-luc constructs. After the synchronization, in the presence of luciferin, light emission was measured. The maximum bioluminescence was set to 0.6. The warm and cold colored lines represent data from two independent biological replicates for the hCry2-luc (−249 to +41) construct and its deletion mutants, respectively. Data sets obtained from two experiments were detrended by subtracting the 24-h running average from raw data.

To investigate whether each of the E-boxes is functional for circadian transcription of the hCry2 gene, we produced constructs carrying hCry2 lacking either one of the E-boxes (Figure 2A). To confirm whether the BMAL1 and CLOCK complex actually activates hCry2 transcription via each E-box, we performed dual luciferase assays with these constructs (Figures 2B). The data showed that overexpression of BMAL1 and CLOCK activated hCry2 transcription via the region spanning from −249 to +41 and that deletion of one of the E-boxes produced severe attenuation of transcription, and deletion of both produced an almost complete loss of activation. Of interest, despite the existence of an intact E1, loss of E2 showed a defect similar to that with double deletion. Although five E-boxes additively act to up-regulate transcription in the case of the regulation of Per1 expression (Hida ), the present data indicate that the tandem E-boxes act to activate hCry2 transcription interdependently. Furthermore, we performed real-time monitoring of cell-autonomous circadian transcription with hCry2-luc constructs lacking either one of the E-boxes (Figures 2C). Of importance, loss of either one led to a complete loss of circadian transcription of the hCry2 gene, suggesting that E1 and E2 play a pivotal role in the endogenous rhythm-generating system in a completely interdependent manner. The single E-box–mediated partial transcriptional activation observed in the reporter assay may be a phenomenon that is detectable only under an overexpression condition. A limitation of these experiments is that, although complete deletion of the E-box (6 base pairs) ensures complete loss of enhancer activity, there may be an increased risk of genomic structural changes, leading to unexpected effects on transcription. The binding assays in Figure 3 were therefore performed using mutant E-boxes with substituted nucleotides rather than deletions.

FIGURE 2:

Interdependent E-boxes enhance nonlinearity in the rhythm-generating system. (A) Schematic representation of deletion mutants of the hCry2 (–249) construct. The deleted sequences are shown at the bottom. (B) Luciferase assay. NIH3T3 cells were transfected with the indicated constructs. Data show relative firefly luciferase activity (Fluc RLU), which was normalized by Renilla luciferase activity (Rluc RLU). Data are represented as mean ± SE for triplicate samples. Student’s t test was performed for statistical analysis (BMAL1-CLOCK [–] versus BMAL1-CLOCK [+]; *p < 0.001, **p < 0.01, ***p < 0.05). (C) U2OS human osteosarcoma cells were transfected with hCry2-luc constructs. After dexamethasone-induced synchronization, in the presence of luciferin, light emission was measured. The maximum bioluminescence was set to 0.6. The warm and cold colored lines represent data from three independent biological replicates for the hCry2-luc (–2925 to +41) construct and deletion mutants of hCry2-luc (–249 to +41), respectively. Data sets obtained from three experiments were detrended by subtracting the 24-h running average from raw data. (D) Transcription function represented by the Hill equation, f(X3, K) = k1/{1 + (X3/K)} (k1 = 1.5, n = 4.2 or 8; dotted line) or its square form, f(X3, K) = [k1/{1 + (X3/K)}]2 (k1 = 1.5, n = 4.2; solid line). Left, CLOCK-BMAL1 level is varied as K ∈ [0, 4], with inhibitor level fixed to X3 = 1. Right, inhibitor level varied as X ∈ [0, 2], with CLOCK-BMAL1 level fixed to K = 1.

FIGURE 3:

The hCry2 tandem E-boxes interact directly with BMAL1 and CLOCK in an interdependent manner. (A) Nucleic acid sequences of a region of hCry2 containing E1 and E2. Asterisks show conserved nucleotides among human, rat, and mouse. The E-boxes are indicated by red squares. Sequences of mutant oligonucleotides used for EMSA and pull-down experiments are shown as mE1, mE2, and mE1 + mE2. A 32P-labeled, double-stranded oligonucleotide that contained the hCry2 tandem E-boxes was used as an EMSA probe. The purified BMAL1 and CLOCK complex was incubated with the probe. The BMAL1 and CLOCK–bound probe is indicated as BM-CL-probe. The complex was shifted by addition of specific antibodies (Super shift). (B) Liver extracts were incubated with double-stranded biotinylated oligonucleotides (hCry2-ODN), which were immobilized on streptavidin-Sepharose beads. Negative control samples were incubated with streptavidin-Sepharose beads without oligonucleotides (Unconjugated). The resulting precipitates were subjected to immunoblot analysis with anti-BMAL1 and anti-CLOCK antibodies. (C) Competitor experiments. The purified BMAL1 and CLOCK complex was incubated with the labeled probe containing the hCry2 tandem E-boxes in the presence of unlabeled probes. The amount of the unlabeled probes is indicated as relative levels to the labeled hCry2 tandem E-box probe. Middle, autoradiography signals quantified with Typhoon FLA9500. Data are representative of more than two independent experiments.

To examine the significance of our experimental findings in silico, we introduced the effect of the two mutually dependent E-boxes to a mathematical modeling. The term k1/{1 + (X3/K)}, known as the Hill equation, describes a sigmoidal binding of transcription factors to a gene promoter. The variable X3 can be considered as nuclear inhibitor concentration (PERs and CRYs in mammals), k1 stands for the production constant, and K represents the level of the BMAL1-CLOCK heterodimer relative to the inhibitor. The steepness of transcription is determined by the Hill coefficient, n. In consideration of our experimental findings, the transcription function in the original equation was modified as the squared form of the Hill equation, [k1/{1 + (X3/K)}]2, so that it represented a multiplicative effect of the transcription through the interdependent binding sites. As shown in Figure 2D, the multiplicative term was able to halve the Hill coefficient that was required to form a sufficient degree of steepness, as well as the nonlinearity of the Hill equation. As shown in both plots, namely, transcriptional activity versus the BMAL1-CLOCK level K (Figure 2D, left) and transcriptional activity versus the nuclear PER-CRY level X (Figure 2D, right), steepness comparable to the original Hill equation with n = 8 was realized by its squared form with only n = 4.2. The original Hill equation with n = 4.2, on the other hand, did not show such sharp steepness. In this way, the multiplicative effect can significantly reduce the Hill coefficient to a realistic value.

Although the reporter assay data show that BMAL1 and CLOCK activate hCry2 transcription via the tandem E-boxes, excluding the possibility that the activation may be due to an artifact indirectly caused by overexpression requires validation of whether circadian transcription factors directly bind to the E-boxes. We therefore performed electrophoretic mobility shift assays (EMSAs) to investigate the physical interaction of purified BMAL1 and CLOCK proteins to double-stranded DNA fragments (Figure 3A). To avoid the possibility that the BMAL1 or CLOCK monomer is structurally unstable and therefore nonfunctional, we coexpressed BMAL1 and CLOCK in HEK293A cells and purified them as a heterodimer. CLOCK expression levels were much lower than those of BMAL1, and a considerable proportion of BMAL1 was therefore purified as a monomer when an antibody against BMAL1 was used. To reduce monomer contamination, we purified BMAL1-saturated CLOCK using an antibody against CLOCK. The purified BMAL1 and CLOCK complex was incubated with a probe containing the hCry2 tandem E-boxes. We confirmed that the electrophoretic mobility of the probe was shifted by the physical interaction with BMAL1–CLOCK. Antibodies against each protein were used to validate that the BMAL1–CLOCK–bound probe was supershifted. A mutation of E1 resulted in a severe attenuation of the affinity of the BMAL1 and CLOCK complex to the probe, and the position of the shifted band was slightly changed, probably because of the reduction of the mole ratio between the binding proteins and probe. A mutation of E2 led to almost complete loss of the physical interaction between them. Taken together, the present EMSA data indicate that each of the hCry2 tandem E-boxes physically interacts directly with the BMAL1 and CLOCK complex in an interdependent and not independent manner. The proximal E-box (E2) may play a role as a scaffold for the distal E-box (E1).

Next we performed hCry2-oligodeoxynucleotide (ODN) pull-down assays and confirmed that endogenous BMAL1 and CLOCK expressed in the liver bound to a fragment containing the hCry2 tandem E-boxes (Figure 3B). The temporal pattern of the binding showed a clear circadian oscillation. The electrophoretically shifted bands corresponded to phosphorylated BMAL1, and both types of BMAL1—phosphorylated and unphosphorylated—bound to the tandem E-boxes. Unexpectedly, a mutation of E1 did not affect the binding affinity of the hCry2-ODN to endogenous BMAL1 and CLOCK, whereas a functional deficiency of E2 caused a severe loss of binding of endogenous circadian transcription factors to the DNA fragment. Unlike the purified protein-based EMSA in Figure 3A, unknown endogenous factors may enhance the affinity of endogenous circadian transcription factors with E2, even in the absence of functional E1.

To further verify the interdependent binding of the hCry2 tandem E-boxes to BMAL1 and CLOCK by another experimental approach, we performed competitor experiments using unlabeled probes that inhibited the recruitment of purified BMAL1 and CLOCK onto the labeled hCry2 tandem E-boxes (Figure 3C). The introduction of the unlabeled tandem E-box probe resulted in a severe reduction of the electrophoretic mobility (lanes 3–5). However, the unlabeled mE1 or mE2 probes had only a slight inhibitory effect on the recruitment of BMAL1 and CLOCK onto the tandem E-boxes (lanes 6–11). In particular, the inhibitory effect of the mE2 probe remained small even when the amount was threefold higher than that of the labeled hCry2 probe. These results confirmed that the loss of either one of the E-boxes led to a severe reduction in binding to BMAL1 and CLOCK. As expected, loss of function of both E-boxes showed no competitor activity even at the high concentration (lanes 12–14). Unlabeled M34 containing a single E-box showed strong competitive activity similar to that of the unlabeled hCry2 tandem E-box probe (lanes 15–17), suggesting that the severe defects observed in mE1 and mE2 competition ability were not simply due to the reduction in the number of functional E-box sites. Although the same mole of unlabeled probes was used for the competition assay, the M34 probe (containing a single E-box) showed slightly stronger competitive activity than the hCry2 tandem E-box probe. This difference in activity may be explained by the fact that M34 was originally identified as a high-affinity DNA-binding site for circadian transcription factors, which was determined by site selection using ∼7 × 107 oligonucleotide randomers (Hogenesch ). However, it is unclear whether the affinity strength of circadian transcription factors to E-box elements is simply proportional to transcriptional activity. Together these results clearly illustrate that functional deficiency of either one of the two E-boxes leads to a severe reduction in the binding affinity of the other E-box with BMAL1 and CLOCK.

To examine the effect of the interdependent binding sites on the autonomous circadian oscillations, we simulated a mathematical model that takes account of the multiplicative effect of the transcription. The simplest, yet-standard mathematical model that captures the essence of transcription-translation feedback loops was developed by Goodwin (1965; Figure 4A, original formula). The variables X1, X2, and X3 can be regarded as mRNA levels of the clock genes, their cytoplasmic protein concentrations, and nuclear inhibitor concentrations, respectively; k1, k2, and k3 represent the production constants; and d1, d2, and d3 determine the degradation rates. The circadian oscillation is generated through transcription function, represented by the Hill equation k1/{1 + (X3/K)}, where the transcription is rhythmically inhibited by the X3 variable through a negative feedback loop. To describe the multiplicative effect of transcription by the CLOCK-BMAL1 dimer through the two E-boxes, we modified the original Goodwin formula as follows (Figure 4A, modified formula):

FIGURE 4:

The multiplicative effect of tandem circadian elements on transcription reduces the Hill coefficient to realistic values. (A) The original and modified Goodwin models. The variables X1, X2, and X3 can be considered as mRNA levels of the clock genes, their cytoplasmic protein concentrations, and nuclear inhibitor concentrations, respectively. k1, k2, and k3 stand for the production constants, and d1, d2, and d3 determine the degradation rates. K represents the level of the CLOCK–BMAL1 complex relative to the inhibitor, and the term 1/{1 + (X3/K)} is known as the Hill equation and describes a sigmoidal binding of transcription factors to gene promoters. The steepness of transcription is determined by the Hill coefficient, n. (B) Left, damped oscillations observed from the original Goodwin model with n = 4.5. From top to bottom, X1, X2, and X3, respectively. The parameter values are set as k1 = 1.5, k2 = 1.5, k3 = 1.5, d1 = 0.15, d2 = 0.15, d3 = 0.15, and K = 1.5. Right, self-sustained oscillations observed from the modified Goodwin model with n = 4.5. From top to bottom, X1, X2, and X3, respectively. The parameter values are set the same as in the original model. (C) Dependence of oscillation amplitude of the X1 variable on the Hill coefficient. Squares and circles correspond, respectively, to the original and modified Goodwin models.

The transcription function in the original equation is thus in a squared form of the Hill equation, such that it represents the multiplicative effect of transcription through the two binding sites. In the original Goodwin model, a Hill coefficient of n > 8 is required to generate a limit-cycle oscillation (Griffith, 1968). However, such an unrealistically large coefficient is often questioned because a Hill coefficient of up to 3 or 4 is usually expected for the formation of protein complexes and their binding to promoters (Gonze and Abou-Jaoude, 2013). In contrast, the modified formula with a squared form of the Hill equation can generate circadian rhythmicity with much smaller Hill coefficients than 8. For instance, Figure 4B shows that only a damped oscillation (no circadian rhythmicity) was observed from the original Goodwin model when the Hill coefficient was set to n = 4.5. In contrast, circadian rhythmicity was easily recovered by the modified model with the same n = 4.5. Figure 4C shows bifurcation diagrams of the original and modified Goodwin models. Dependence of the oscillation amplitude of the X1 variable on the Hill coefficient is drawn. Here the oscillation onset point is known as the Hopf bifurcation point. Although the original Goodwin model shows oscillation onset at n = 8, the modified Goodwin model gives rise to a self-sustained oscillation at n = 4.2. This is mathematically reasonable because a squared form of the Hill equation can enhance nonlinearity of the transcription term, as shown in Figure 2D, which supports oscillations. The multiplicative term can halve the Hill coefficient required to form a sufficient degree of nonlinearity of the Hill equation. In this way, the multiplicative effect can significantly reduce the Hill coefficient to realistic values. We further carried out a bifurcation analysis by computing the largest eigenvalue of the Jacobian matrix of our model and confirmed that the oscillation onset is due to supercritical Hopf bifurcation. This bifurcation type was also identified in an experimental observation of oscillation onset in cyanobacteria circadian clocks (Murayama ).

The multiplicative effect of transcription can be further expanded to RORE, a pair of which is known to activate and suppress transcription of Bmal1 in a mutually dependent manner via the action of ROR and REV-ERB, respectively (Preitner ; Akashi and Takumi, 2005). We therefore extended the mathematical model to take into account the function of the two interdependent ROREs. Because the binding of ROR to RORE is mutually exclusive with that of REV-ERB, we omitted introduction of the ROR action to simplify the model (Figure 5A). Following Becker-Weimann , we added a BMAL1 feedback loop to the core loop of PER-CRY in the mathematical model (Figure 5B). Here the BMAL1-CLOCK het­erodimer activates the transcription of Per, Cry, and Rev-erb genes. After translation and posttranslational processes, PERs and CRYs inhibit their own transcription, forming the main negative feedback loop. In the BMAL1 loop, transcription of Bmal1 leads to the formation of nuclear BMAL1, which, as a complex with CLOCK, restarts the transcription of Per, Cry, and Rev-erb genes. Because REV-ERB represses the transcription of Bmal1, BMAL1 indirectly inhibits its own transcription, forming an additional negative feedback loop. As implemented in the Goodwin model in Figure 4, the effect of transcription through interdependent binding sites was described in multiplicative terms for both E-boxes and ROREs (Figure 5A). As shown in Figure 5C, the multiplicative terms effectively reduced the minimum Hill coefficient needed to sustain the circadian oscillation. Whereas the minimum Hill coefficient was n = 5.3 in the case in which no multiplicative effects were considered, it was reduced to n = 2.6 when the multiplicative effect was considered for E-boxes and further reduced to n = 2.55 when the multiplicative effects were considered for both E-boxes and ROREs. The oscillation amplitude expanded rapidly even though the minimum Hill coefficient was not substantially reduced. The contribution of interdependent ROREs to the enhancement of circadian amplitude was thus confirmed.

FIGURE 5:

The multiplicative effect of tandem ROREs further reduces the Hill coefficient. (A) Equations for the interconnected feedback model (top) and the modified model (bottom). X1, X2, and X3 represent concentrations of Per and Cry mRNA, PER/CRY complex in the cytoplasm, and PER/CRY complex in the nucleus, respectively. X4, X5, and X6 represent concentrations of Bmal1 mRNA, cytoplasmic BMAL1 protein, and BMAL1 protein in the nucleus, respectively. X7 and X8 stand for concentrations of BMAL1-CLOCK heterodimer and REV-ERB protein, respectively. The parameter values were set as v1b =18 nM h–1, k1b = 1 nM, K1 = 0.65 nM, c = 0.001 nM, s = 5, k1d = 0.12 h–1, k2b = 0.3 nM–1 h–1, q = 2, k2d = 0.07 h–1, k2t = 0.24 h–1, k3t = 0.02 h–1, k3d = 0.12 h–1, v4b = 3.0 nM h–1, k4b = 1 nM, K4 = 0.9 nM, u = 0.5, k4d = 1.8 h–1, k5b = 0.14 h–1, k5d = 0.03 h–1, k5t = 0.15 h–1, k6t = 0.06 h–1, k6d = 0.03 h–1, k6a = 0.03 h–1, k7a = 0.003 h–1, k7d = 0.02 h–1, v8b = 10.6 nM h–1, k8b = 1 nM, K = 1.1 nM, r = 1, v = 2, and k8d = 1.5 h–1. (B) Schematic illustration of the mathematical model. BMAL1–CLOCK heterodimer activates the transcription of Per, Cry, and Rev-erb genes. After translation and posttranslational processes, nuclear PER/CRY complex inhibits transcription of Per, Cry, and Rev-erb genes, forming the main negative feedback loop. Through REV-ERB, which represses the transcription of Bmal1, BMAL1 indirectly inhibits its own transcription, forming an additional negative feedback loop. (C) Dependence of the Hill coefficient and the oscillation amplitude of X1 variable on the interdependence of E-boxes and ROREs. A multiplicative effect was considered for transcription of only E-boxes (purple), only RORα (green), both E-boxes and RORα (blue), or none (red).

DISCUSSION

Many previous studies identified components of the mammalian circadian clock and revealed the framework of its molecular mechanism. Nevertheless, gaining a comprehensive understanding of the system for generating circadian oscillation by biological experimental approaches alone is difficult. To aid understanding, experimental data from biological studies have been frequently introduced into mathematical models. This approach has uncovered several issues that were not raised by the fragmentary information obtained by biological experiments. The Goodwin model is widely used to simulate mammalian circadian transcription-translation feedback loops in silico. However, to realize the degree of nonlinearity required for autonomous transcriptional oscillation, the Hill coefficient must be adjusted upward to experimentally unrealistic values. Although this issue can be mathematically avoided by assuming high nonlinearity in several circadian processes, evidence for the presence of sufficient nonlinearity in these processes has not been reported. We therefore conducted the present study to identify other biological processes that contribute to the enhancement of nonlinearity.

Although recent comprehensive experimental approaches, such as DNA microarray, chromatin immunoprecipitation sequencing (ChIP-seq), and RNA sequencing, have substantially aided our understanding of the framework of the mammalian circadian system, detailed validation of the molecular mechanism of circadian transcription remains to be done in some clock genes, including Cry2, a paralogue of Cry1, which encodes one of the core components of the mammalian clock. In vivo studies have demonstrated that both CRY1 and CRY2 play essential roles in circadian period control, with their deficiencies resulting in different effects on period length (van der Horst ). In contrast, ex vivo and in vitro studies have revealed that CRY1, but not CRY2, is indispensable for maintaining circadian amplitude (Liu ; Westermark ). Three sets of evidence may explain the different phenotypes observed between Cry1 and Cry2 deficiency. First, the circadian phase of Cry1 expression is remarkably delayed compared with that of other circadian transcriptional repressor genes (Per1, Per2, Per3, and Cry2; Matsumura ). Although this might underlie the different roles of the two CRYs in circadian period regulation, it is unclear whether it can explain why Cry2 is dispensable for circadian amplitude. Second, intracellular protein levels of Cry1 are nearly double those of Cry2 (Lee ). Although this suggests that the contribution of Cry2 protein to the circadian oscillator is smaller than that of Cry1, it might not explain the normal circadian amplitude observed in Cry2-deficient cells. Finally, in vitro reporter assays have revealed that Cry1-mediated transcriptional inhibition is much stronger than that by Cry2 (Akashi ), suggesting that Cry1 could strongly compensate for a Cry2 deficiency. This is likely the most convincing reason for the dispensability of Cry2 in circadian amplitude regulation.

Traditional and classical molecular biological and biochemical analyses of circadian transcription might provide insights into biological processes that enhance nonlinearity in the circadian system. Comprehensive ChIP-seq data have already indicated the physical interaction of the 5′-flanking region of the hCry2 gene with the circadian transcription factors BMAL1 and CLOCK (Koike ), but the functional significance of the interaction remains undefined. In the present study, we adopted an experimental approach that enables identification of the promoter regions required for circadian expression based on actual oscillatory phenotypes and identified a 17–base pair sequence required for rhythmic expression. Furthermore, our data indicate that this small region indeed contains two E-boxes, which are indispensable not only for the physical interaction with the circadian transcription factors BMAL1 and CLOCK but also for actual circadian transcription of the hCry2 gene.

The enhancer element E-box is not specific to clock gene transcription, raising the issue of what defines the functional difference between circadian E-boxes and those controlling the expression of genes related to other diverse biological processes. Although studies on the binding specificity of transcription factors to the E-box on target alleles have provided evidence for the contribution of sequences surrounding E-boxes (Munoz ; Wang ), the mechanism determining the functional characteristic of circadian E-boxes remains undefined. The present study indicates that interdependent transcriptional regulation via tandem E-box elements may be a hallmark characteristic of circadian oscillation in transcription. A tandem E-box motif was also found in the 5′-flanking region of the Per and Dbp genes (Akashi ; Nakahata ; Rey ), suggesting that interdependent transcriptional regulation may be a common mechanism in clock gene expression. Of importance, we previously reported that two ROREs interdependently control circadian expression of the Bmal1 gene: loss of the distal RORE function leads to almost complete functional loss of proximal RORE (Akashi and Takumi, 2005). These two ROREs are tandemly repeated and spaced by 25 nucleotides (Preitner ). These previous studies indicated that the role of these tandem motifs is to enhance circadian amplitude. In the present study, however, we speculated that the interdependent behavior of these motifs might contribute not only to the amplitude but also to the limit-cycle oscillation in clock gene expression. According to our mathematical model, interdependent tandem E-boxes and ROREs can indeed enhance nonlinearity in the circadian transcription of clock genes. Interdependent transcriptional regulation via tandem elements may therefore provide sufficient nonlinearity not only to the PER/CRY main loop but also to the BMAL1 subloop in the circadian feedback system, indicating that this nonlinear regulation may be a common and hallmark characteristic of transcription of core clock genes.

Our interpretation of the functional role of interdependent tandem E-boxes and ROREs is as follows. In the case that two elements act independently and additively, the plot of transcription versus transcriptional activator concentration shows low nonlinearity. In contrast, when two elements act interdependently, the plot shows a steep sigmoidal pattern because the activation of transcription requires the binding of transcriptional activators to both. Transcriptional efficiency of interdependent elements should be lower than that of independent elements during the circadian phase when the nuclear amount of transcriptional activators is small but should increase in a steep sigmoidal manner as transcriptional activators accumulate in the nucleus. We mathematically confirmed this scenario by modifying the Hill equation to a squared form. An increase in nonlinearity was detected in the plot of transcription versus activator concentrations K1 and K4 (BMAL1 and CLOCK for E-box, ROR for RORE), and this ultrasensitivity was similarly observed in the plot of transcription versus inhibitor concentrations X3 and X8 (PER and CRY for E-box, REV-ERB for RORE). Taken together, interdependent transcriptional regulation via tandem elements may play an indispensable role in realizing cosine-like expression rhythms of clock genes. Furthermore, interdependent tandem enhancers on multiple clock genes might cooperatively enhance nonlinearity in the whole circadian feedback system, which would lead to limit-cycle oscillations in clock gene expression.

It is important to note that the role of the cooperative effect of multiple factors on transcriptional regulation has been expressed mathematically (as quadratic, cubic, and higher-order terms) in other contexts, including the general theory of gene expression and regulation (Bintu ), immunological memory imprinting in T-helper lymphocytes (Hofer ), delay equation modeling of circadian regulatory networks (Korencic ), and the 12-h rhythmicity observed in various gene expressions (Westermark and Herzel, 2013). In gene regulation, a number of interacting genes, proteins, and molecules often coexist to perform a particular function. The cooperative effect might be a common biological strategy to enhance overall functionality and potentially overcome weak regulatory functions of individual genes.

MATERIALS AND METHODS

Plasmid construction

A BAC clone (RP11-618K13) containing the transcriptional regulatory region for the human Cry2 gene was purified to generate an hCry2 promoter–driven luciferase expression vector. The hCry2 region from –2925 to +41 (+1 is a putative transcription start site) was amplified by PCR using primers containing restriction enzyme sites that were required for subsequent subcloning. The primer sequences were as follows: hCry2-BamHI, 5′ GCA GGA GGA TCC CTT GAG 3′; and hCry2-NcoI, 5′ CCA TGG CTG TCC AGA CTG CTC CAG 3′. The PCR product was TA-cloned into the pCRII TOPO vector (Thermo Fisher Scientific), and the insert sequence was verified. Subsequently the insert was excised using BamHI and NcoI and subcloned into the pGL3-basic vector (Promega).

A deletion series of the transcriptional regulatory region for hCry2 was constructed by combining double digestion and self-ligation procedures. The restriction enzymes used for double digestion of hCry2(-2925∼+41)-pGL3 were as follows: hCry2(-1904), KpnI (two restriction sites, one in the vector multicloning site and one in the insert); hCry2(-1014), XhoI and PstI; hCry2(-500), XhoI and SnaBI; and hCry2(-249), XhoI and SacII.

Deletion mutants of hCry2(-249)-pGL3 were generated by inverse PCR using back-to-back primers to amplify the entire plasmid vector, and linear PCR products were subsequently self-ligated to form circular DNA. The entire sequence of the transcriptional regulatory region was verified in each deletion construct to confirm the presence of the designed deletion and the absence of unexpected PCR errors.

The coding regions for mouse BMAL1 and CLOCK were amplified by nested PCR using mouse liver cDNA. After TA-cloning and sequence verification, the coding regions were excised using restriction enzymes in the vector multicloning site and subcloned into the pcDNA3 expression vector (Invitrogen), as described in a previous study (Akashi ).

Cell culture and transfection

Cells were grown in DMEM (supplemented with antibiotics and 10% fetal bovine serum [FBS]) and cultured in 5% CO2. Transfection of plasmid DNA into cells was performed with FuGENE 6 (Promega) or Lipofectamine plus (Invitrogen).

Real-time monitoring of luciferase activity in living cells

U2OS human osteosarcoma cells were cultured and transfected with indicated combinations of expression vectors and incubated for 20 h. Two hours after the treatment with dexamethasone (50 nM), the medium was replaced with fresh culture medium (supplemented with 1% FBS). In the presence of 0.1 mM luciferin, light emission was measured and integrated for 1 min at intervals of 15 min with an LM2400 photomultiplier tube (Hamamatsu Photonics). The bioluminescence data were processed through the 24-h moving average, and the maximum level was set to 0.6.

Luciferase reporter assay

NIH3T3 cells were cultured in 24-well plates and transfected with indicated combination of vectors (a promoter-driven firefly luciferase vector, a Renilla luciferase vector as an internal control, plus clock protein expression vectors). About 20 h after transfection, cell lysates were prepared and used according to the manufacturer’s protocol of the Dual Luciferase Assay System (Promega). Luminescence was measured with GloMax 20/20 luminometer (Promega). Data (Fluc/Rluc relative light units [RLUs]) show relative firefly luciferase activity (Fluc RLU), which was normalized by the Renilla luciferase activity (Rluc RLU). The Student’s t test was used for comparisons between two groups.

Preparation for EMSA

Proteins were purified from HEK293A cells transfected with expression plasmids. Immunoprecipitation was performed with antibody-conjugated beads, and the immunoprecipitates were eluted with an excess amount of tag peptide. Purification was confirmed by staining SDS–PAGE gels. A double-stranded oligonucleotide was used as an EMSA probe. A 25-pmol amount of the double-stranded oligonucleotide was labeled with [γ-32P]ATP using a polynucleotide kinase (Toyobo). G-50 Spin Oligo Columns (GE Healthcare) were used to remove free ATP from labeled probes.

EMSA

Binding buffer contains 15 mM Tris-HCl (pH 7.5), 15 mM NaCl, 1.5 mM EDTA, 1.5 mM dithiothreitol, 7.5% glycerol, 0.3% NP-40, and 1 μg/ml bovine serum albumin. A radiolabeled probe was added at 0.5 nM, and purified clock proteins were mixed and incubated for 10 min. Finally, antibodies recognizing the clock proteins were added to perform supershift experiments. All of the reactions were incubated at 25°C. The reactions were loaded onto 4% nondenaturing polyacrylamide gels and resolved at 150 V for 2 h at 4°C. The gels were dried, and autoradiography signals were quantified with Typhoon FLA9500 (GE Healthcare Life Sciences). The nucleotide sequence of the M34 competitor was GGGACACGTGACCC (Hogenesch ).

Animals

Animals were maintained on a 12:12-h light–dark cycle (light on at 9:00 A.M.), and allowed ad libitum access to food and water. Male C57Bl/6J mice were purchased from Japan SLC (Shizuoka, Japan). All experiments were in accordance with the rules of the Yamaguchi University Animal Usage Committee.

Pull-down experiment

Liver extracts were prepared by homogenizing in ice-cold incubation buffer. Then the extracts were incubated with a double-stranded biotinylated oligonucleotide that had been immobilized on streptavidin-Sepharose beads (Amersham). After washing with the incubation buffer, the resulting bound protein was subjected to immunoblot analysis with anti-BMAL1 and anti-CLOCK antibodies.