Authors: Lanxiao Zhu, Bin Qian, Xin Zhang, Lin-Lin Hu
Categories: Original Research, circadian rhythm, cardiovascular disease, heart, multi-omics analysis, translation efficiency, BMAL1/REV-ERBα
Source: Vascular Health and Risk Management
Doi: 10.2147/VHRM.S578582
Authors: Lanxiao Zhu, Bin Qian, Xin Zhang, Lin-Lin Hu
Cardiovascular diseases remain the leading cause of death and economic burden worldwide. Increasing evidence indicates that sleep disturbance and circadian rhythm disruption are major risk drivers for hypertension, coronary artery disease, heart failure, and arrhythmia. Although the classical transcription–translation feedback loop (TTFL) model explains the basic mechanism of rhythm generation, increasing evidence suggests that the heart—an organ with high metabolic demand—maintains circadian stability through coordinated transcriptional, translational, and post-translational regulation.
We developed an integrative, time-resolved, multilayer in silico framework to systematically analyze cardiac circadian regulation by combining mouse heart time-series RNA-seq (GSE54650), proteomics (PXD002870), phosphoproteomics (PXD036824), BMAL1 and Rev-erbα ChIP-seq, and enhancer RNA (eRNA) datasets. Rhythmicity was assessed using MetaCycle, with cross-layer comparisons evaluating concordance and divergence between transcriptomic and proteomic rhythms, and translation efficiency (TE) estimated from protein-to-mRNA ratios. Enhancer–gene coupling, transcription factor binding, and phosphorylation motif analyses were integrated to investigate multilayer regulatory coordination.
We identified 2552 rhythmic transcripts and 139 rhythmic proteins, with only 31 genes rhythmic at both layers, indicating substantial RNA–protein phase decoupling in the heart. Temporal stability of TE correlated positively with protein amplitude, suggesting that stable translation supports robust protein rhythmicity. Phosphoproteomic analyses revealed enrichment of SP motifs mediated by proline-directed kinases in rhythmic proteins. BMAL1 binding was associated with enhanced transcriptional amplitude, whereas REV-ERBα binding was associated with delayed target gene expression, forming complementary enhancer-level regulatory dynamics.
This study supports a multilayered integrative model of cardiac circadian regulation in which rhythmic gene expression is jointly shaped by transcriptional activation, translational precision, and post-translational modification. By extending the classical “clock–transcription–protein” paradigm, our findings highlight enhancer-level control mediated by BMAL1 and Rev-erbα as an important mechanism contributing to the stabilization of cardiac circadian timing.
Cardiovascular diseases remain the leading cause of death and economic burden worldwide, making their prevention and control a major challenge for public health systems.1,2 Growing evidence indicates that sleep disorders and circadian rhythm disruption are important drivers of cardiovascular risk, although their underlying mechanisms remain incompletely understood.3 Despite substantial advances in chronobiology, how the heart—a highly metabolic organ—maintains rhythmic stability through coordinated regulation across multiple molecular layers remains unresolved. Early studies established that circadian rhythms are driven by a transcription–translation feedback loop (TTFL) involving core factors such as BMAL1, CLOCK, PER, CRY, and REV-ERBα,4–6 laying the molecular foundation of the field. With the development of multi-omics and systems biology, researchers have increasingly recognized that cardiac rhythmicity is shaped not only by transcriptional control but also by translational regulation, post-translational modifications, and metabolic feedback as coupled, multilayer mechanisms.3,7,8 Elucidating the molecular basis of cardiac circadian rhythms may help explain how insomnia and related circadian disturbances act as independent cardiovascular risk factors, and may provide a conceptual framework for the application of chronotherapy in cardiovascular disease prevention and treatment.3,9
Although multi-omics studies have revealed layer-specific molecular regulatory programs during cardiac development, the temporal coordination of these layers remains insufficiently defined.10 Evidence across tissues indicates that phase lags and amplitude mismatches between RNA and protein are common, suggesting that cross-layer rhythms may coexist in coupled and desynchronized states.8,11 At the chromatin level, rhythmic BMAL1 binding is associated with transcriptional amplitude, whereas REV-ERBα is often linked to phase-delayed transcriptional repression. Beyond transcription, post-translational modifications such as phosphorylation are thought to influence the stability of rhythmic proteins and circadian signal propagation.12–17 In highly dynamic tissues such as the heart, protein stability and turnover are increasingly recognized as additional contributors to circadian output beyond transcription and translation alone.18 However, how these regulatory layers coordinate over time to jointly determine the phase and magnitude of rhythmic gene expression remains insufficiently understood from a cross-layer, systems-integrative perspective.
Among core circadian regulators, BMAL1 and REV-ERBα were selected in this study not simply because of their canonical roles within the TTFL, but because they represent two functionally complementary modes of circadian regulation.17 BMAL1 primarily acts as a transcriptional activator that enhances rhythmic amplitude,16,19,20 whereas REV-ERBα mediates phase-delayed transcriptional repression.21 Importantly, both factors bind enhancer regions and exert direct regulatory control over circadian timing.16,19,20 This activation–repression pairing provides a mechanistically informative framework to investigate how circadian signals are propagated, temporally gated, and reshaped across transcriptional, translational, and post-translational layers in the metabolically demanding cardiac environment.
Prior studies at the molecular, systems, and clinical levels have elucidated the circadian features and regulatory mechanisms of the cardiovascular system. Basic research shows that multiple functions of the heart and vasculature (eg, blood pressure, heart rate, myocardial contractility, and energy metabolism) exhibit circadian oscillations driven jointly by the central suprachiasmatic nucleus (SCN) and peripheral clocks.6 Circadian disruption and sleep disorders are associated with substantially increased risks of hypertension, coronary heart disease, heart failure, and arrhythmias,3,15 whereas maintaining normal circadian rhythms helps preserve the cardiovascular system’s “physiological resilience”.3 Extensive experimental evidence indicates that impaired BMAL1 function is linked to abnormalities such as disrupted myocardial calcium homeostasis, heightened oxidative stress, and apoptosis, thereby promoting pathologies including diabetic cardiomyopathy and myocardial infarction.22,23 REV-ERBα deficiency abolishes the circadian rhythm of fatty-acid oxidation and causes metabolic mismatch, leading to myocardial energy crisis and heart-failure phenotypes.21 Moreover, multi-omics studies reveal that circadian regulation extends well beyond transcription, encompassing layers such as translational efficiency, phosphorylation, and metabolic rhythms,8,10,24 and forming phase-coordinated rhythmic networks across organs.8,25 Collectively, these studies support a regulatory continuum linking the molecular clock, multilayer metabolism, organ function, while highlighting the need for integrative temporal analyses.
Despite substantial progress, important gaps remain. Most studies focus on a single molecular layer or specific disease models, lacking a cross-layer, time-series, systems-integrative perspective. Although the TTFL model explains the generation of rhythms in principle, its coupling with post-translational oscillators (PTOs) and metabolic feedback has not been fully elucidated.26 In addition, the temporal coordination mechanisms and phase differences across layers (transcription, translation, phosphorylation) lack quantitative characterization, and applications of multi-omics integrative analysis are still at an early stage.8,10,24 This emphasis on integrative and stratified frameworks parallels contemporary cardiovascular risk models, such as the PRECISE-HBR score for post-PCI bleeding risk prediction, in which multiple biological or clinical dimensions are combined to improve interpretability and decision-making.27 Although chronotherapy has been proposed as a potential intervention, the translational pathway from molecular rhythms to clinical practice remains to be further explored.9
To address these gaps, we performed a time-resolved, systems-level integration of multi-omics datasets from the mouse heart, including RNA-seq, proteomics, phosphoproteomics, BMAL1 and REV-ERBα ChIP-seq, and enhancer RNA (eRNA) profiles. Using BMAL1 and REV-ERBα as regulatory anchors, we integrated these datasets along the temporal axis to elucidate the cross-layer coupling among transcriptional, translational, and post-translational regulatory tiers; clarify the complementary regulatory roles of BMAL1 and REV-ERBα at circadian enhancers; and propose a multilayer analytical framework for cardiac circadian regulation, providing new avenues to understand the molecular basis of rhythmic homeostasis in the heart and its potential interventional relevance in cardiovascular diseases associated with circadian disruption.
Detailed information on all transcriptomic, proteomic, phosphoproteomic, and regulatory datasets used in this study is summarized in Table 1. Table 1Overview of Multi-Omics Datasets Used in This StudyProfileLibraryTime PointReferenceAnimal InformationMature RNARNA-seq24 (2 h intervals over 48 h)GSE54650 (GEO)12Mouse; C57BL/6; male; 6 weeks; n=3Translation rateEstimated (RNA/protein)24 (inferred values)Computed from matched datasets/ProteomicsMass spectrometry (MS)24 (2 h intervals over 48 h)PXD002870 (PRIDE @ ProteomeXchange)28Mouse; multiple strains; male; 9–12 weeks of age; 2 mice per experimental time point in the original turnover studyPhosphoproteomicsTMT10-plex1 time pointPXD036824 (PRIDE @ ProteomeXchange)10Not reportedEnhancer RNAGRO-seq1 time pointENCFF665YOG29Not reportedBmal1ChIp-seq1 time point (ZT6)GSM3003981 (GEO)20Mouse; C57BL/6J; male; 3–6 months old; n = 3Rev-erbChIp-seq2 time point (ZT4/ZT16)GSE153150 (GEO)21Not reported
RNA expression data were obtained from the GSE54650 dataset,12 which comprises 24 time points from mouse heart tissue sampled every 2 hours over a 48-hour period (ZT18–ZT64). Probe-level measurements were extracted from the original expression matrix (series.txt), and probe IDs were mapped to gene symbols using the platform annotation file. For genes with multiple probes, expression values were averaged. All gene expression time series were standardized using Z score normalization and subsequently used as input for circadian rhythm analysis.
Protein expression data were obtained from the PXD002870 dataset,28 which includes mouse heart samples collected every 2 hours over a 48-hour period (ZT0–ZT46; 24 time points). Proteins with more than 7 missing values or zero expression at more than half of the time points were excluded. The remaining protein expression profiles were normalized using Z-score transformation and subsequently subjected to circadian rhythm analysis.
To assess the relationship between transcriptional and protein-level rhythmicity, proteins were categorized into three groups on the basis of their overlap with rhythmic 11 Both – rhythmic at both the RNA and protein levels;Gain – rhythmic only at the protein level;Lost – rhythmic only at the RNA level.
Functional enrichment analysis was performed independently for each group.
To investigate the role of phosphorylation in the maintenance, acquisition, and loss of circadian rhythmicity, we first performed a joint classification of genes on the basis of time series transcriptomic and proteomic data.11 Genes were categorized into three Both (rhythmic at both the RNA and protein levels), Gain (rhythmic only at the protein level), and Loss (rhythmic only at the RNA level).
Phosphoproteomics data were obtained from the PXD036824 dataset in the ProteomeXchange repository, which is derived from 8-week-old (P56) mouse heart tissue.10 Given that these data represent a single time point, we focused on phosphorylation intensity and motif enrichment across the three gene groups to infer the regulatory contributions of posttranslational modifications during rhythmic state transitions.10,30 This approach aimed to elucidate how phosphorylation contributes to both the stability and plasticity of the circadian regulatory system.31
BMAL1 ChIP-seq data from mouse heart tissue at Zeitgeber time 6 (ZT6) were obtained from the GEO dataset GSM3003981 (GEO).20 The peak file (GSM3003981_Heart_BMAL1_ChIP_All_Replicates.bed.gz) was annotated using the ChIPseeker R package32 with the reference transcript database TxDb.Mmusculus.UCSC.mm10.knownGene. Each peak was assigned to the nearest transcription start site (TSS) within ±3 kb to define putative target genes.
To identify BMAL1-bound enhancers, ChIP-seq peaks were overlapped with previously defined rhythmic enhancer regions. Enhancers were classified as BMAL1_bound = TRUE or FALSE on the basis of whether they overlapped a ChIP peak. Given that BMAL1 binding typically leads to transcriptional activation with a 2–6 hour delay,16 we defined a “daytime enhancer window” as rhythmic enhancers with peak phases between ZT8 and ZT14.32,33 This window was used to identify BMAL1-bound enhancers that are likely functionally active.
These enhancers were intersected with a set of 31 dual-layer rhythmic genes identified from integrated RNA and protein rhythmicity analysis to define BMAL1-regulated rhythmic genes.11 To quantify BMAL1 binding strength, the corresponding BigWig signal file was used to calculate the mean signal intensity within ±3 kb of each peak centre. Spearman correlation analysis was then performed between BMAL1 binding intensity and the amplitude (AMP) of associated rhythmic genes.
Putative enhancer regions were defined using H3K27ac ChIP-seq peak data from adult mouse heart tissue provided by the ENCODE project (accession: ENCFF665YOG)29 The peak file was imported into R as a GRanges object and annotated by computing the nearest distance to transcription start sites (TSSs) using the TxDb.Mmusculus.UCSC.mm10.knownGene reference transcript database.
Enhancer–gene pairs were established on the basis of proximity to the TSS. If the paired gene was classified as rhythmic by MetaCycle,34 the corresponding enhancer was designated a rhythmic enhancer. For each enhancer region, the average phastCons score (60-way vertebrate conservation)35 was calculated to assess evolutionary conservation.
Enhancers with a phastCons_score ≥ 0.7 and target genes peaking in the ZT8–ZT14 window (phase24 ∈ ZT8–ZT14)20 were defined as daytime-typical rhythmic enhancers. Genes associated with these enhancers were subjected to GO and KEGG pathway enrichment analysis using clusterProfiler v4.0.36 Enhancers linked to nonrhythmic genes were used as a control group. To compare conservation levels between rhythmic and nonrhythmic enhancers, a Wilcoxon rank-sum test was performed on the distribution of phastCons scores.
BigWig files of Rev-erbα at ZT4 and ZT16 were downloaded from GSE153150 (GEO).21 Using rtracklayer37 and GenomicRanges,38 the average signal of each rhythmic enhancer was extracted at both time points, and ΔBinding = ZT4 − ZT16. Enhancers with ΔBinding > 0.0189 (top 5% of the overall difference distribution) are defined as Rev-bound enhancers, while the rest are classified as unbound. We subsequently compared the two groups of enhancers in terms of AMP, phastCons scores, and phase differences using the Wilcoxon rank-sum test to assess significance and visualize the results with ggplot2.
The rhythmicity of RNA expression and protein abundance was evaluated using the meta2d algorithm in the MetaCycle R package,34 with the period range set to 20–28 hours. For RNA-level analysis, genes whose adjusted p value was < 0.05 were considered significantly rhythmic. For protein-level analysis, owing to the lower detection power and signal-to-noise ratio, a raw p value < 0.1 was used as the criterion to identify candidate rhythmic proteins. For estimated TE rhythmicity, an adjusted p value < 0.1 was applied.
Using proteomic data (PXD002870, 48 hours with sampling every 2 hours) and matched mRNA expression data (GSE54650), the translation efficiency (TE) was calculated as TE = log2(protein/mRNA).39,40 Only genes with detectable expression in both datasets and no missing (NA) or infinite values were retained for analysis. For each gene, we calculated the mean TE (TE_mean), standard deviation (TE_sd), and coefficient of variation (TE_CV = TE_sd/TE_mean). To mitigate the influence of low expression levels or noise, genes whose TE_mean was in the bottom 10% of the overall distribution or whose TE_CV was less than 0.3 were filtered out. As no ribosome profiling (Ribo-seq) data were available for the analyzed time series, TE values derived here represent an estimated proxy rather than a direct measurement of translational activity. Accordingly, analyses involving TE focus on relative temporal stability and variability rather than absolute translational rates, and related interpretations should be made in this context.39,40
Welch’s t test was used to compare TE_mean and TE_CV across different gene categories, with statistical significance defined as p < 0.05.
Phase information (peak expression time in ZT) of rhythmic genes was extracted using the MetaCycle package.34 Bar plots and polar coordinate diagrams were generated to visualize the phase distributions of RNA and protein rhythmicity separately, allowing the observation of clustering patterns across the circadian cycle.
For genes exhibiting rhythmicity at both the RNA and protein levels, the phase difference (Phase_Diff = Protein_Phase − RNA_Phase) was calculated. Histograms were constructed to evaluate the distribution of phase differences, indicating overall synchrony, delays, or advances in protein expression relative to RNA expression.
Genes were further grouped on the basis of the direction of the phase difference (delay or advance), and functional enrichment analysis was performed using clusterProfiler v4.036 to identify Gene Ontology (GO) biological processes associated with each group.
Finally, Spearman correlation analysis was conducted between the RNA and protein phases to assess the temporal coupling between transcription and translation.
To comprehensively analyse circadian gene expression in mouse cardiac tissue, we analysed RNA-seq data from the GSE5465012 dataset and performed rhythmicity analysis using the MetaCycle (meta2d) algorithm.34 The cycle range was set between 20 and 28 hours, with an adjusted p value (adj.p.) threshold of less than 0.05 for significance. A total of 2552 genes exhibiting significant circadian rhythmicity were identified (Figure 1A). Figure 1(A) Heatmap of the top 50 high-amplitude rhythmic genes in the mouse heart. RNA-seq expression values were Z-score normalized across circadian time points and visualized as a heatmap. Warmer colors indicate higher relative expression and cooler colors indicate lower expression. Genes were hierarchically clustered to highlight temporal expression patterns and phase relationships. (B) Peak phase distribution of rhythmic transcripts in the mouse heart. A polar histogram showing the distribution of rhythmic genes (n = 2552) across circadian phases (0–24 h) at the RNA level. Bar height represents the number of genes peaking at each phase.
A polar clock plot (Figure 1B) revealed that the peak expression phases of circadian genes were unevenly distributed throughout the 24-hour cycle, indicating distinct rhythmic clustering. Notably, phase 18 (corresponding to 00, or early evening) was the dominant peak, with approximately 200 genes reaching their maximum expression at this time. Secondary peaks were observed at phase 21 and phase 15, whereas relatively few genes peaked at Phase 3 or Phase 9, indicating a nonuniform phase distribution across the 24-hour cycle.
To further elucidate the biological roles of the 2552 circadian genes, we performed GO and KEGG enrichment analyses using clusterProfiler v4.0.36
The GO enrichment results revealed that these circadian genes were significantly associated with biological processes (BP), including muscle cell differentiation, muscle cell proliferation, response to hypoxia, cardiac muscle tissue development, and regulation of vasculature and angiogenesis (Figure 2A). KEGG pathway analysis revealed that rhythmic genes are enriched in several key signalling pathways, including the circadian rhythm pathway, MAPK and PI3K-Akt signalling pathways, HIF-1 signalling, fluid shear stress and atherosclerosis, and metabolic and endocrine-related pathways, such as glucagon signalling and the AGE-RAGE signalling pathway, in diabetic complications (Figure 2B).
Figure 2Functional enrichment analysis of circadian transcripts. (A) Gene Ontology (GO) analysis of 2552 rhythmic genes showing enrichment in biological processes related to muscle cell differentiation, hypoxia response, and vascular development. The top 10 terms are shown based on −log2(p.adjust) values. (B) KEGG pathway analysis showing enrichment of rhythmic transcripts in signalling and metabolic pathways, including circadian rhythm, MAPK, PI3K–Akt, HIF-1, and glucagon signalling. Bars represent the top 10 enriched pathways ranked by −log2(p.adjust).
To assess whether transcriptional rhythmicity is propagated to the protein level, we analysed mouse cardiac proteomic data from the PXD002870 dataset,28 which includes samples collected every 2 hours over a continuous 48-hour period. Circadian rhythm detection was performed using the meta2d algorithm in the MetaCycle R package,34 with the period range set to 20–28 hours. Owing to the relatively low signal-to-noise ratio at the proteomic level, rhythmic proteins were identified using a raw meta2d p value threshold of < 0.1. This analysis revealed 139 proteins that exhibit significant circadian rhythmicity. After annotation to gene symbols, these rhythmic proteins corresponded to 279 unique genes, a markedly smaller number than that observed at the transcriptomic level.
Among genes with valid measurements at both the RNA and protein levels, only 31 genes exhibited significant circadian rhythmicity at both layers, accounting for approximately 0.63% of the total analyzed genes. This limited overlap indicates that rhythmic transcription alone is insufficient to ensure rhythmic protein expression in the mouse heart.
To compare the temporal distribution of rhythmic signals across molecular layers, we generated a radar plot dividing the 24-hour cycle into six consecutive Zeitgeber time (ZT) windows (ZT0–4, ZT4–8, ZT8–12, ZT12–16, ZT16–20, and ZT20–24) (Figure 3A). Rhythmic transcripts displayed a relatively uniform phase distribution with a prominent peak at ZT12–16. In contrast, rhythmic proteins were enriched in the ZT4–8 and ZT12–16 intervals and depleted in the ZT8–12 and ZT16–20 intervals, showing clear differences in phase distributions between RNA and protein layers. Figure 3Phase characteristics and cross-layer comparisons of rhythmic genes at the RNA and protein levels. (A) Radar plot showing the distribution of rhythmic genes across six circadian phase windows at the RNA (Orange) and protein (blue) levels. Values represent the number of rhythmic genes in each Zeitgeber time window. (B) Histogram of phase differences between protein and RNA peak phases (ΔPhase = protein phase − RNA phase) for 31 dual-rhythmic genes. Positive values indicate protein-lagging genes (red), whereas negative values indicate RNA-lagging genes (blue).
To further characterize cross-layer temporal relationships, we analyzed phase differences for the 31 genes rhythmic at both levels. The phase difference (ΔPhase = Phase_protein − Phase_RNA) was calculated and adjusted to the interval [–12, +12]. RNA and protein phases exhibited only a weak and nonsignificant correlation (Spearman’s ρ = 0.10, p > 0.05).
Analysis of ΔPhase distributions revealed that 51.61% (n = 16) of these genes exhibited protein expression peaks lagging behind RNA peaks by approximately 5–8 hours (Figure 3B). Functional enrichment analysis of these lagging genes did not identify significant GO terms or KEGG pathways.
Conversely, 48.39% (n = 15) of the genes exhibited earlier protein peaks relative to RNA peaks (ΔPhase < 0). Functional enrichment analysis revealed significant enrichment of metabolic processes, including lipid biosynthesis, gluconeogenesis, hexose biosynthesis, and acyl-CoA metabolism (Supplementary Figure S1A). In contrast, the lagging group showed limited enrichment, with only the spliceosome pathway being significantly overrepresented (Supplementary Figure S1B).
Together, these analyses revealed extensive RNA–protein phase differences in the mouse heart, indicating limited concordance between transcriptional and proteomic circadian rhythms.
To further investigate whether protein-level rhythmicity depends on mRNA rhythmicity, we selected genes whose expression was valid at both the RNA and protein levels (n = 4958) and classified them into four groups on the basis of a previously described cross-layer rhythmicity classification 11 Rhythmic-Both: significant rhythmicity at both the RNA and protein levels (n = 31, ~0.63%);Rhythm-Lost: significant rhythmicity at the RNA level but not at the protein level (n = 928, 18.72%);Rhythm-Gained: no significant rhythmicity at the RNA level but significant at the protein level (n = 109, 2.20%);Nonrhythmic: no significant rhythmicity at either level (n = 3890, 78.45%).
This distribution showed that a large proportion of genes were classified as Rhythm-Lost, whereas only a small subset retained rhythmicity at both RNA and protein levels (Figure 4). Based on this classification, subsequent analyses focused on the three rhythm-related groups. Figure 4Distribution of rhythmic gene classes across RNA and protein layers. Bar plot showing the classification of 4958 genes with detectable expression at both RNA and protein levels into three rhythm-related Rhythmic-Both (n = 31), genes rhythmic at both RNA and protein levels; Rhythm-Lost (n = 928), genes rhythmic only at the RNA level; and Rhythm-Gained (n = 109), genes rhythmic only at the protein level. Genes nonrhythmic at both layers (~78%) are not shown.
GO enrichment analysis revealed that Rhythmic-Both genes were significantly enriched in a single GO term, regulation of lipid biosynthetic process (Supplementary Figure S2A). In contrast, Rhythm-Lost genes were enriched in processes such as the regulation of protein-containing complex assembly, the regulation of mRNA metabolic processes, protein folding, the regulation of translation, and muscle cell differentiation (Supplementary Figure S2B). Rhythm-Gained genes were enriched in GO terms related to protein targeting, cytoplasmic translation, protein–RNA complex assembly, protein localization, and fatty acid oxidation (Supplementary Figure S2C). KEGG pathway enrichment analysis revealed that Rhythmic-Both genes were significantly enriched in the spliceosome pathway (Supplementary Figure S3A). Rhythm-lost genes were significantly enriched in pathways related to cellular structure, stress response, and protein processing, including the cytoskeleton in muscle cells, fluid shear stress and atherosclerosis, PI3K-Akt signalling, and protein processing in the endoplasmic reticulum (Supplementary Figure S3B). Rhythm-Gained genes did not show significant KEGG enrichment, likely reflecting their smaller number and functional heterogeneity.
Together, these results demonstrate that cross-layer circadian rhythmicity in the heart follows structured and non-random patterns, underscoring the role of posttranscriptional and posttranslational regulation in shaping protein rhythmicity.
To explore the role of translation efficiency (TE) in the formation of protein rhythmicity, we estimated the TE by calculating the log2 ratio of protein to RNA expression levels at 24 time points per gene as a proxy for translation efficiency. Genes were classified into three Rhythmic-Both, Rhythm-Gained, and Rhythm-Lost.
No significant difference was observed in the mean TE (TE_mean) across the three groups, suggesting that overall translation efficiency is not the primary determinant of protein rhythmicity.
In contrast to TE_mean, temporal fluctuations in TE, measured by the TE standard deviation (TE_sd, Kruskal–Wallis p = 6.93×10^−5^) and the coefficient of variation (CV_TE, p = 0.00235), differed significantly across groups. The Rhythm-Lost group presented the highest TE_sd and CV_TE, corresponding to the largest temporal variability in translation efficiency. Conversely, the Rhythmic-Both group displayed the lowest TE_sd and CV_TE, while the Rhythm-Gained group had intermediate values.
Further pairwise comparisons revealed that the Rhythm-Lost group had significantly greater CV_TE than both the Rhythmic-Both and Rhythm-Gained groups did. In contrast, genes in the Rhythmic-Both and Rhythm-Gained groups exhibited much lower CV_TE values (Figure 5). Figure 5Translation efficiency (TE) dynamics across rhythmic gene groups. Violin and box plots showing the distributions of (A) mean translation efficiency (TE_mean), (B) standard deviation of TE (TE_sd), and (C) coefficient of variation of TE (CV_TE) in Rhythmic-Both, Rhythm-Lost, and Rhythm-Gained gene groups. TE_mean did not differ significantly among groups (Kruskal–Wallis test, p = 0.625), whereas TE_sd and CV_TE showed significant differences (p = 6.93×10^−5^ and p = 0.00235, respectively).
To further assess the relationship between translation efficiency variability and protein rhythmicity, Spearman correlations were calculated between CV_TE and protein expression amplitude (AMP). Negative correlations were observed in all three groups. The strongest inverse association was observed in the Rhythm-Gained group (ρ = −0.575, p = 4.13×10^−11^), followed by the Rhythmic-Both group (ρ = −0.566, p = 0.0011), whereas the Rhythm-Lost group showed a moderate but significant negative correlation (ρ = −0.334, p < 2.2×10^−16^) (Figure 6). Figure 6Association between translation efficiency variability and protein rhythmic amplitude. Scatter plots show CV_TE versus protein expression amplitude (AMP) for the Rhythmic-Both (left), Rhythm-Gained (middle), and Rhythm-Lost (right) groups. Spearman correlation analysis showed significant negative correlations in all three groups. The red line indicates the linear regression fit; Spearman ρ and p values were calculated separately.
To investigate the role of phosphorylation in circadian rhythm modulation, we classified phosphosites into three groups based on the rhythmic status of their corresponding Rhythmic-Both, Rhythm-Gained, and Rhythm-Lost (n = 92, n = 267, and n = 3296 phosphosites, respectively). The Rhythm-Lost group contained the largest number of sites, whereas the Rhythmic-Both group contained relatively few sites. Analysis of normalized phosphorylation intensities showed higher phosphosite levels in the Rhythm-Gained and Rhythm-Lost groups compared with the Rhythmic-Both group, while the corresponding protein abundance in the Rhythm-Lost group was lower (Wilcoxon rank-sum tests; pairwise p values = 0.063–0.37; Figure 7A).
Figure 7Phosphorylation-level differences and motif signatures across protein rhythmicity classes. (A) Violin and box plots showing normalized phosphosite intensity (log2-transformed phosphosite intensity relative to mean protein abundance) in the Rhythm-Gained, Rhythm-Lost, and Rhythmic-Both groups. (B–D) Heatmaps showing amino-acid frequency distributions in 13-mer sequence windows centered on phosphosites (positions −6 to +6) for the Rhythm-Gained (B), Rhythm-Lost (C), and Rhythmic-Both (D) groups. (E) Sequence logos of the same 13-mer windows for the Rhythm-Gained, Rhythm-Lost, and Rhythmic-Both groups, highlighting differences in local motif composition around phosphorylation sites.
To characterize sequence features underlying differential phosphorylation patterns, motif analysis was performed by constructing 13-mer amino acid sequences centred on each phosphosite using sequence windows provided by MaxQuant.41 Visualization revealed that the Rhythm-Gained and Rhythmic-Both groups exhibited clear enrichment of proline (P) at the +1 position (the 8th amino acid), forming a canonical SP motif, whereas the Rhythm-Lost group showed a more dispersed amino acid distribution at this position, resulting in less defined motif boundaries (Figure 7B–D). Fisher’s exact test further confirmed that the SP motif was significantly enriched in the Rhythmic-Both (p = 0.08985) and Rhythm-Gained (p = 0.0051) groups compared with the Rhythm-Lost group (Figure 7E).
To quantify the relationship between BMAL1 binding and transcriptional amplitude of circadian genes, ChIP-seq signal intensity at ZT6 was averaged over promoter regions (TSS ± 3 kb) and mapped to RNA rhythmic genes identified by MetaCycle. BMAL1 promoter occupancy showed a weak but statistically significant correlation with RNA expression amplitude (ρ = 0.079; p = 0.00016; Figure 8A). To evaluate BMAL1 effects at the enhancer level, we compared AMP distributions between genes associated with BMAL1-bound versus nonbound daytime enhancers (ZT8–ZT14). Genes linked to BMAL1-bound enhancers exhibited significantly higher AMP values (Wilcoxon rank-sum test, p = 8.37 × 10^−21^; Figure 8B). In addition, BMAL1 binding strength showed a mild association with RNA peak phase (ρ = 0.046; p = 0.029; Supplementary Figure S4). Figure 8Association between BMAL1 binding and transcriptional amplitude of circadian genes. (A) Spearman correlation between BMAL1 ChIP-seq signal intensity at gene promoters (TSS ± 3 kb) and RNA expression amplitude (AMP) for circadian genes. A weak but significant positive correlation was observed (ρ = 0.079, p = 0.00016). The red line indicates the linear regression fit; Spearman ρ and p values were calculated separately. (B) Violin and box plots comparing RNA expression amplitude (AMP) between genes associated with BMAL1-bound enhancers and non-BMAL1-bound enhancers during the daytime. Genes linked to BMAL1-bound enhancers exhibited significantly higher AMP values (Wilcoxon rank-sum test, p = 8.37×10^−21^), suggesting that enhancer-level BMAL1 binding amplifies rhythmic transcriptional output.
To further examine how BMAL1-associated transcriptional regulation propagates to the protein level, we screened for genes harboring BMAL1-bound enhancers within the defined diurnal enhancer window and intersected these genes with those exhibiting rhythmic expression at both the mature RNA and protein levels. This approach identified two representative targets, Mcf2l and Nfia, which were subjected to detailed analysis.
Time-series expression analysis revealed that Mcf2l exhibited a mature RNA peak at approximately ZT12, whereas its protein abundance peaked at ZT20–ZT22, corresponding to a delay of approximately 8–10 hours (Figure 9A). In contrast, Nfia showed nearly synchronous RNA and protein expression peaks (Figure 9B). Figure 9BMAL1-centred regulatory cascade and posttranscriptional delay analysis of Mcf2l and Nfia. (A) Time-series RNA and protein expression profiles of Mcf2l, showing a pronounced delay (~8–10 h) in protein accumulation relative to the RNA peak. (B) Time-series RNA and protein expression profiles of Nfia, with closely aligned phases, consistent with predominantly transcriptional regulation. (C) Comparison of total RNA-to-protein phase delay between Mcf2l and Nfia. (D) Decomposition of phase delay into RNA→TE (Orange) and TE→protein (blue) components, indicating a translation efficiency–mediated delay in Mcf2l but not in Nfia. (E) Schematic timeline of Mcf2l regulation, illustrating BMAL1 binding at ZT6, enhancer activation at ZT8–10, RNA peak at ZT12–14, TE increase at ZT16–18, and protein peak at ZT20.
Although multitimepoint ribosome profiling data were not available, translational efficiency (TE) was estimated using the protein-to-RNA ratio. The inferred TE trajectory of Mcf2l exhibited a delayed peak relative to its RNA expression (Figure 9C and D), whereas Nfia displayed a relatively stable TE profile, consistent with the close alignment of its RNA and protein phases. A schematic timeline summarizing the temporal relationships among RNA expression, protein abundance, and inferred TE for Mcf2l is shown in Figure 9E.
Based on H3K27ac ChIP-seq data for mouse heart tissue provided by the ENCODE project, more than 100,000 active enhancer regions were identified and linked to their nearest gene transcription start sites (TSSs) to generate an enhancer–gene pair dataset. By intersecting these enhancer-associated genes with RNA-level circadian genes identified by MetaCycle, a total of 51,834 enhancers linked to rhythmic genes were obtained and defined as candidate rhythmic enhancer regions.
To assess evolutionary conservation, phastCons60way cross-species conservation scores were applied. Rhythmic enhancers exhibited a broad range of conservation values (Supplementary Figure S5), and statistical comparison revealed significantly greater conservation in rhythmic enhancers than in nonrhythmic enhancers (Wilcoxon rank-sum test, p = 7.3×10^−8^).
Further filtering identified enhancer regions with conservation scores ≥ 0.7 whose target genes peaked between ZT8 and ZT14 as high-confidence rhythmic enhancer candidates, yielding 205 enhancer–gene pairs. KEGG pathway enrichment analysis showed that genes associated with these enhancers were significantly enriched in the circadian rhythm pathway (Supplementary Figure S6A). GO Biological Process enrichment analysis further identified terms related to circadian regulation of gene expression, cardiac muscle cell contraction, and potassium ion transmembrane transport (Supplementary Figure S6B).
Among typical rhythmic enhancers whose linked target genes peaked during the daytime window (ZT8–ZT14; n = 16,640), 0.82% (n = 137) overlapped BMAL1 ChIP-seq signals (Supplementary Figure S7).
To investigate the regulatory role of REV-ERBα in circadian enhancers within the mouse heart, we analysed ChIP-seq data collected at ZT4 and ZT16 and calculated the differential binding score for each enhancer (ΔBinding = ZT4 − ZT16). Following a dynamic ChIP-seq time-point comparison strategy inspired by Koike et al (2012), enhancers with ΔBinding > 0.0189 (top 5%) were defined as REV-ERBα-bound.
Compared with genes associated with unbound enhancers, genes associated with REV-ERBα-bound enhancers exhibited significantly greater transcriptional amplitude (AMP) (p < 2.2e−16, Wilcoxon test, Figure 10A), indicating an association between REV-ERBα binding and transcriptionally active circadian regulatory elements. Additionally, these enhancers exhibited significantly higher sequence conservation (phastCons score) (p = 1.96e−10; Figure 10B). Figure 10Comparison of Rev-erbα-bound and unbound enhancer regions. (A) Violin plot comparing RNA expression amplitude (AMP) of genes associated with Rev-erbα-bound versus unbound enhancers; Rev-erbα-associated genes exhibit higher transcriptional amplitudes. (B) Violin plot showing the distribution of phastCons conservation scores for Rev-erbα-bound and unbound enhancers, indicating higher evolutionary conservation of Rev-erbα-bound regions. (C) Violin plot comparing circadian phase distributions of genes linked to Rev-erbα-bound and unbound enhancers, showing that Rev-erbα-associated genes tend to peak later in the circadian cycle.
To further assess the temporal characteristics associated with REV-ERBα binding, we compared the phase distributions of target gene expression. Phase analysis of target gene expression revealed that the expression peaks corresponding to Rev-bound enhancers were generally delayed, occurring mainly between ZT12 and ZT20, whereas the expression peaks of unbound enhancers were concentrated between ZT6 and ZT12 (p < 2.2e−16; Figure 10C).
When the relationship between the estimated translation efficiency (TE) and Rev-erbα binding was examined, no significant difference in the estimated TE distribution was observed between the two groups (p = 0.7611, Wilcoxon test).
This study demonstrates that circadian regulation in the mouse heart is characterized by widespread RNA–protein phase decoupling rather than a simple linear propagation of rhythmic signals from transcription to protein expression. Although rhythmic transcription is prevalent, only a limited subset of genes retains rhythmicity at the protein level, and RNA–protein phase concordance is generally weak, indicating that such decoupling represents a fundamental regulatory feature rather than technical noise. Similar RNA–protein uncoupling has been reported across multiple tissues, where rhythmic protein expression frequently exhibits both leading and lagging patterns relative to transcript oscillations, reflecting multilayered posttranscriptional and posttranslational regulation rather than transcriptional control alone.24,42 Together, these findings suggest that cardiac protein rhythmicity is actively reshaped downstream of transcription to fine-tune protein activity in accordance with circadian-driven physiological demands.7 Conceptually, this multilayered regulatory architecture aligns with integrative frameworks commonly applied in cardiovascular research to improve system-level interpretabilityg.27
The prevalence of rhythm-lost and rhythm-gained gene classes indicates a fundamentally nonlinear architecture linking transcriptional and proteomic circadian regulation in the heart. The frequent loss of protein rhythmicity despite oscillatory transcription suggests active posttranscriptional filtering rather than regulatory failure, underscoring the limited predictive power of transcriptional oscillations alone.11,24 Conversely, rhythm-gained genes demonstrate that rhythmic protein expression can arise independently of RNA rhythmicity, likely through translational and posttranslational control mechanisms.31,42 These findings support a selective routing model in which transcriptional rhythms are differentially transmitted, reshaped, or newly generated at downstream regulatory layers rather than uniformly propagated from RNA to protein.7 From a functional perspective, genes that lost protein rhythmicity were preferentially associated with processes related to cardiac metabolism, protein homeostasis, and structural organization. These functional categories are closely linked to the high and temporally regulated energy demands of the heart,6,43 suggesting that posttranscriptional filtering may selectively disrupt rhythmic protein output in pathways that are particularly sensitive to translational control and proteostatic balance.24,42 In contrast, rhythm-gained genes were more frequently linked to translational and protein complex–related processes, consistent with the notion that rhythmic protein expression can be actively generated downstream of transcription.31,42 Collectively, these results indicate that nonlinear circadian signal transmission across molecular layers is functionally biased toward specific biological processes relevant to cardiac physiology.
Our results indicate that the temporal stability of translation efficiency, rather than its mean level, is a key determinant of protein rhythmicity. Excessive temporal fluctuations in translation efficiency impair the generation or maintenance of rhythmic protein expression, whereas stable translation supports robust protein oscillations. This interpretation is consistent with prior proteomic studies showing that translational variability and noise can disrupt circadian protein accumulation independently of transcriptional amplitude.42,44 Together, these findings suggest that translational noise, rather than insufficient translational capacity, underlies the loss of protein rhythmicity in many transcriptionally rhythmic genes.
Importantly, protein rhythmicity is not determined by translational regulation alone but reflects the balance between protein synthesis and degradation. In highly metabolically active tissues such as the heart, rapid protein turnover does not necessarily lead to unstable protein levels; instead, fast synthesis–degradation cycles may coexist with buffering mechanisms that preserve temporal protein profiles. Although direct measurements of protein half-life were not available in this study, prior work indicates that protein stability and degradation kinetics can substantially modulate circadian protein oscillations independently of transcriptional rhythms.42,44 In this context, excessive temporal variability in translation efficiency may amplify degradation-driven noise, whereas stable translation may buffer against rapid turnover and support sustained protein rhythmicity.
Beyond translational regulation, posttranslational modification further contributes to the stabilization of rhythmic protein expression. Rhythmic proteins exhibited enrichment of the canonical SP motif, implicating proline-directed kinases in sustaining or establishing protein rhythmicity, consistent with prior reports of circadian phosphorylation dynamics in clock-regulated proteins.30,31 Such conserved kinase recognition motifs may facilitate coordinated regulation of protein activity and stability, providing an additional buffering layer that enables rhythmic protein expression even in the absence of rhythmic transcription.42
Enhancer-level analyses indicate that BMAL1 and REV-ERBα exert complementary but distinct roles in shaping transcriptional rhythms in the heart. BMAL1 primarily amplifies transcriptional amplitude through enhancer engagement,16,20 whereas REV-ERBα delays transcriptional timing by temporally repressing enhancer activity.21,45 However, neither factor alone dictates protein rhythmicity. Representative targets illustrate that rhythmic transcription established at enhancers can be differentially transmitted to the proteome depending on downstream translational and posttranslational regulation.24,42 Together, these findings support a hierarchical model in which enhancer-level clock control establishes a permissive transcriptional framework that is subsequently filtered by downstream regulatory layers.7
This study has several limitations. First, translation efficiency was estimated indirectly from matched RNA and protein abundance due to the lack of ribosome profiling data, which may not fully capture instantaneous translational dynamics. Second, protein turnover rates were not directly measured. In addition, the phosphoproteomic data used in this study were derived from single-time-point measurements, which limits direct assessment of the temporal dynamics of phosphorylation rhythms. Given the high metabolic activity of the heart, protein rhythmicity likely reflects the integrated effects of translational regulation and protein degradation. Future studies incorporating ribosome profiling, time-resolved phosphoproteomics, and direct measurements of protein turnover will be necessary to further refine the multilayer regulatory framework proposed here. Finally, as this study was performed in wild-type mice, it remains to be determined how disease-associated mutations, including those affecting sarcomeric proteins, may influence multilayered circadian regulation in pathological cardiac settings.
In conclusion, this study demonstrates that cardiac circadian gene expression is regulated in a multilayered and asynchronous manner, rather than through a simple linear propagation from transcription to protein abundance. Distinct subsets of genes achieve rhythmic protein expression through different regulatory routes operating across transcriptional, translational, and post-translational layers.
At the transcriptional level, BMAL1-dependent enhancer activation primarily amplifies rhythmic mRNA expression during the daytime, but the transmission of these rhythms to the proteome is frequently delayed and shaped by downstream translational control. In parallel, REV-ERBα-mediated repression acts in advance at circadian enhancers to impose temporal delays on transcriptional activation, thereby fine-tuning the phase of rhythmic gene expression. Beyond transcription, temporal fluctuations in translation efficiency emerge as a critical determinant of whether rhythmic mRNA output is converted into stable protein oscillations. Finally, phosphorylation-dependent post-translational regulation, particularly involving proline-directed kinases, appears to support the stabilization and maintenance of rhythmic protein expression.
Together, these findings refine the classical clock–transcription–protein framework by emphasizing the essential contributions of post-transcriptional and post-translational regulation in shaping cardiac circadian output. This multilayered regulatory model provides a conceptual foundation for understanding rhythm-associated cardiovascular dysfunctions and may inform future investigations into rhythm-associated cardiovascular dysfunctions and chronotherapy-based cardiovascular interventions.6,43