Multi-omics profiling reveals rhythmic liver function shaped by meal timing

The multi-omics study robustly identified tens of thousands of PTM sites and proteins. These included 10,874 proteins, 30,463 phospho-sites, 11,136 ubiquityl-sites, 7944 succinyl-sites, and 5245 N-glycosyl-sites, of which 17,359 phospho-sites, 8430 ubiquityl-sites, 7446 succinyl-sites, and 4305 glycosyl-sites are quantified (<30% missing values across samples) (Fig. 1b). In addition, we quantified 571 lipid species covering 29 lipid classes in the lipidomics study, and 50 lipid species covering 4 lipid classes that belong to N-Acyl-Phosphatidylethanolamine (NAPE), endocannabinoids (EC) and acyl amino acids (acyl-AA) (Fig. 1b).

We sought to identify PTM sites on the circadian clock proteins. In total, 57 sites on circadian clock proteins that can be modified by phosphorylation or ubiquitylation were identified in this multi-omics study (Supplementary Table 1), including the recently reported mouse PER2 Ser693/Ser697 and Cryptochrome 1 (CRY1) Ser588 sites²². For example, Rev-Erbα/NR1D1 lysine 456 was identified as a ubiquitylation site (Supplementary Fig. 1a). Lysine 456 mediates the physical interaction between NR1D1 and PER2, the loss of which reduces the interaction by 60% (ref. ⁴²). Rev-Erbα/NR1D1 Ser279 was identified as a potential phosphorylation site among three candidate sites (Supplementary Fig. 1b), the mutation of which reduces the interaction with the FBXW7 E3 ubiquitin ligase and subsequent degradation⁴³. Notably, PER2-pSer971 had the most peptide spectrum matches (Supplementary Fig. 1c), but the in vivo regulation of this modification has not been elucidated.

To determine daily rhythms of proteins and lipids in response to TRF, we applied a stringent algorithm. The multi-omics profiles were subjected to three different rhythmicity detection methods, including a nonparametric detection method (Rhythmicity Analysis Incorporating Nonparametric method, or RAIN), a parametric method (CircaCompare), and a hybrid method (MetaCycle)^44–46. Diurnal rhythmicity is determined by adjusted P-values < 0.05 in RAIN and MetaCycle and a P-value < 0.05 in Circacompare. Consequently, we detected robust diurnal rhythms in 45% (2139/4778) of phospho-proteins, 38% (483/1277) of ubiquityl-proteins, 34% (1983/5773) of proteins, 16% (196/1252) of N-glycosyl proteins, and 1.9% (24/1286) of succinyl-proteins in the mouse livers (Fig. 1b). The phospho-proteome appears much more rhythmic than the other forms of PTMs, which is consistent with the high degree of circadian rhythms of the phospho-proteome in the liver of ad libitum-fed mice²². In contrast, protein succinylation exhibits the lowest degree of diurnal rhythm, at least within the depth of the current proteome profiling.

Considering the great percentage of diurnal proteins in the phospho-proteome, we examined diurnal rhythms in the mouse phospho-proteome first. Most diurnal phospho-proteins oscillated around ZT4 under DRF and at ZT8 and ZT20 under NRF (Supplementary Fig.). The amplitude and relative amplitude (the ratio between amplitude and baseline) of diurnal phospho-proteins were significantly lower under DRF compared to NRF (Supplementary Fig.). Pathway enrichment based on the phase of diurnal phospho-proteins revealed that almost all rhythmic pathways peak in the sleep phase (ZT0-8) during DRF and concentrate in the first half of the daytime (ZT3-6) (Supplementary Fig.). This is represented by the regulation of protein targeting (ZT3), cell matrix adhesion (ZT4), response to calcium ion (ZT5) and macroautophagy (ZT6) (Supplementary Fig.). In the NRF group, enriched pathways almost exclusively peaked in the second half of the sleep phase (ZT6-11) (Supplementary Fig.) and included lymphocyte homeostasis (ZT7), regulation of cell matrix adhesion (ZT8), mitochondrial membrane organization (ZT9) and drug catabolic process (ZT10) (Supplementary Fig.). 2a 2a 2b 2c 2b 2d

We performed pathway analysis based on phase-shifts and found that most rhythmic pathways were shifted in phase for 4–5 h by DRF, compared to NRF (Fig. 2d). Rhythmic pathways that were poorly entrained by meal time included inositol lipid-mediated signaling (phase-shift 3 h, Kuiper test, q-value = 0) and receptor-mediated endocytosis (phase-shift 4 h, Kuiper test, q-value = 0) (Fig. 2e). We found that the receptor-ligand machinery involved in lipid clearance exhibited robust diurnal rhythms under both DRF and NRF, but it was not altered in phase between DRF and NRF, as indicated by diurnal phosphorylation profiles of LRP1 (LDL receptor-related protein 1), LDLRAP1 (LDL receptor adaptor protein 1) and LRPAP1 (LRP associated protein 1) (Fig. 2f).

The most food-entrainable pathways exhibited a phase-shift of 6–9 h and included monocarboxylic acid biosynthetic process (phase-shift 9 h), regulation of circadian rhythm, fatty acid metabolism, mitochondrial matrix and transcriptional corepressor activity (Fig. 2g). Notably, phosphorylated forms of the core clock proteins, such as PER2 and Rev-Erbα/NR1D1, were inverted in phase by DRF compared to NRF (Fig. 2h). PER2-pSer971 exhibited robust diurnal rhythms in the livers of both DRF and NRF mice, and its diurnal abundance was completely reversed in phase by DRF compared to NRF. Overall, 23.1% (143/620) of the dual-cycling phospho-proteins were entrained by DRF in mouse livers compared to NRF, despite a robust complete entrainment of the circadian clock by meal timing.

Next, we explored diurnal rhythms of the hepatic proteome. Proteomics by a data-independent acquisition mode provided us with a deep, robust, and quantitative approach to analyze diurnal rhythms among 5773 hepatic proteins. Most diurnal proteins peaked around ZT20 and ZT0 during DRF and NRF, respectively (Supplementary Fig. 3a). The relative amplitude is slightly smaller under DRF compared to NRF (Wilcoxon test P = 0.026), but the amplitude is not statistically different between DRF and NRF (Supplementary Fig. 3a). Notably, 28.0% or 380 of the 1358 diurnal proteins from DRF mice peaked at ZT20. In line with this, enriched DRF diurnal pathways peaked around ZT20 (Supplementary Fig. 3b), including xenobiotic metabolic process (ZT18), serine family amino acid metabolism (ZT19), cellular hormone metabolic process (ZT19), neurotransmitter metabolic process (ZT20) and monosaccharide biosynthetic/catabolic process (ZT21) (Supplementary Fig. 3c). Enriched NRF diurnal pathways peaked from ZT23 to ZT1 (Supplementary Fig. 3b), including xenobiotic metabolic process (ZT22), regulation of cholesterol biosynthetic process (ZT23), Protein N-linked glycosylation (ZT23), lipid storage (ZT24 or ZT0) and peroxisome organization (ZT1) (Supplementary Fig. 3d). These features seen in livers from NRF mice are found in ad libitum-fed mouse livers collected in constant darkness, as well²¹.

We examined diurnal rhythmicity of representative proteins in the lipid metabolic pathway (Circacompare method, P < 0.05), and found that the rhythmicity of the lipid droplet coat protein Perilipin-2 (PLIN2) was shifted in phase for 7.5 h and increased significantly in baseline (mesor) by DRF compared to NRF (Fig. 3f). DRF abolished the rhythmicity of PLIN3 protein, whereas it shifted the phase of PLIN5 protein by 6.8 h without altering the baseline (mesor) and amplitude (Fig. 3f). Proteins involved in fatty acid desaturation (FADS2), very-long/long chain fatty acyl coenzyme A biosynthesis (SLC27A2, ACSL4) and oxidation (ACOX2) and degradation of bioactive lipids (FAAH) were all increased in baseline abundance by DRF in the liver and shifted in phase for various hours (Fig. 3f). Overall, lipid metabolic proteins exhibit a pathway-level increase in abundance and a phase-shift of 8 h.

Diurnal phase analysis of the di-glycine peptide-enriched ubiquityl-proteome profiles performed in our study revealed that diurnal ubiquityl-proteins from DRF mice were evenly distributed throughout the diurnal cycle, whereas diurnal ubiquityl-proteins from NRF were mainly in the dark/active phase in mouse livers (Supplementary Fig. 4a). Both the average relative amplitude and the average amplitude were decreased by DRF compared to NRF (Supplementary Fig. 4a). The 26S proteasome activity was diurnal in livers from DRF mice and peaked around ZT20.5 (ZT19.7 from Circacompare method, adjusted P = 0.0076; ZT21 from RAIN P = 0.022); however, it was not rhythmic in livers from NRF mice (Supplementary Fig. 4b). This is associated with a significant decrease in the baseline (mesor) activity.

Pathway enrichment analysis revealed that rhythmic pathways under DRF reached peak time hours from ZT10 to ZT20, spanning from the late feeding/day phase to the fasting/night phase (Supplementary Fig.). Diurnal signatures found in the hepatic ubiquityl-proteome from DRF mice are represented by lipid/steroid metabolic process (ZT10), xenobiotic metabolic process (ZT11), carbohydrate metabolic process (ZT13) and response to lipid (ZT14) (Supplementary Fig.). In contrast, a diverse array of pathways is rhythmic and peak around ZT17 in the ubiquityl-proteome from NRF mice (Supplementary Fig.), including alpha/cellular amino acid metabolic process (ZT16), lipid homeostasis (ZT17), lipid modification (ZT18) and response to metal ion (ZT19) (Supplementary Fig.). 4c 4d 4c 4d

To identify diurnal features driven by DRF, we took a supervised N-integration sparse discriminant analysis and integrated diurnal multi-omics data from the same mouse liver samples. Proteins, phospho-proteins, and ubiquityl-proteins with daily rhythms in at least one TRF condition served as the input data to mixOmics analysis^47,48. The modeling applies the N-integrative supervised analysis with DIABLO, because the same N samples are measured on different ‘omics platforms. We divided the input data into three subsets based on the time of day; i.e., morning subset (sampled at ZT0 and ZT4, n = 16 per group), evening subset (sampled at ZT8 and ZT12, n = 16 per group) and night subset (sampled at ZT16 and ZT20, n = 16 per group), and fed into the statistical modeling. The results showed that the evening subset returned the best-fit model with an overall balanced error rate (BER) of 0.0312 compared to the other two subsets (morning subset BER = 0.278 and night subset BER = 0.203).

Next, we focused on the mixOmics model based on the evening subset of multi-omics data and found that the statistical model could clearly segregate TRF groups in the phospho-proteome and proteome dataset, and to a lesser degree in the ubiquityl-proteome (Fig. 4f). Importance plotting shows the rank and weight of the diurnal proteins that distinguish DRF from NRF. We found that phosphorylated heme-binding protein 1 (p-HEBP1) and phosphorylated Rev-Erbα/NR1D1 (p-NR1D1), which is a nuclear receptor for heme⁴⁹, and p-PER2 contributed the most (Fig. 4g).

A list of 20 diurnal proteins from the proteomics dataset contributes to the ability to distinguish the effects of DRF from NRF on the liver diurnal rhythm (Fig. 4g). Enzymes involved in fatty acid metabolism are also represented, such as ACSL4, ACOX2, and SLC27A2. Strikingly, four of the five contributing proteins in the diurnal ubiquityl-proteome (CYP2C29, SLC27A2, CYP2C50, and CYP4A14) are involved in arachidonic acid metabolism, which is part of fatty acid metabolism (Fig. 4g). To reveal the correlation among the contributing features, we generated a circular correlation plot and found that phospho-proteins and ubiquityl-proteins are tightly correlated (Fig. 4h). In particular, p-PER2 is positively correlated with ubiquityl (ub)-SLC27A2 whereas the latter is negatively correlated to p-NR1D1.

We further confirmed this feature by examining the clustered image map, where p-PER2 and many of the fatty acid metabolic enzymes, e.g., ACSL4, SLC27A2 and its ubiquityl form, ubiquityl forms of CYP2C29, CYP2C50, and CYP4A14, show greater weight in DRF compared to NRF (Fig. 4i). Overall, our integrative model indicates that the circadian clock, particularly represented by p-PER2, and its tightly associated proteins in the fatty acid metabolism pathway, is a key diurnal feature in response to TRF in the mouse liver.

Considering that a rather small percentage of succinyl-proteins cycled in the livers of TRF mice (Fig. 1b), we examined diurnal rhythms among individual sites of succinylation. We found that diurnal succinyl-sites are clustered into two groups that peak around ZT10 and ZT20, respectively, under DRF, whereas succinyl-sites are generally scattered around the clock except in one peak at ZT16 under NRF (Supplementary Fig. 5c). There is no significant difference between DRF and NRF, with respect to the average amplitude (Supplementary Fig. 5d). In total, there are 97 diurnal succinyl-sites, seven sites of which are shared between DRF and NRF mice (Fig. 5g). Consistent with the fact that protein succinylation occurs mainly in mitochondria⁵⁰, almost all of the dual-cycling succinyl-sites were found in mitochondrial proteins, except K414 of CYP2D26 (Fig. 5h). Rhythmic succinylation at K644 of mitochondrial trifunctional enzyme alpha subunit HADHA and K93 of SUCLA2 were inverted in phase by DRF compared to NRF, and these are the only two food-entrainable succinyl-sites identified in this study (Fig. 5h, i).

The phase-inverted lipids consist mainly of Bis-monoacylglycero-phosphate (BMP) and glycerophospholipids (Fig. 6b). BMP is primarily present in endosomal and lysosomal compartments⁵¹. BMP species peaked in the sleep and wake time under DRF and NRF, respectively, in synchrony with the feeding time (Fig. 6c). Acylcarnitine species were either moderately phase-shifted or phase-locked in the liver from DRF mice compared to NRF (Fig. 6b, c). Remarkably, DRF induced robust diurnal rhythms among more than 80 species of triacylglycerols (TAGs) and 8 species of diacylglycerols (DAGs) in the liver (Fig. 6d and Supplementary Fig. 6). In contrast, NRF entrained diurnal rhythms of one TAG species in the liver, while it entrained robust diurnal rhythms of acylcarnitines - essential for fatty acid oxidation - and phosphoglycerolipids, primarily phosphatidylserines (PSs) (Fig. 6d).

We matched the profiles of TAGs and DAGs with those of monoacylglycerols (MAGs) and biophysically relevant lipids by targeted quantification profiling. Diurnal rhythms in N-acylphosphatidyl ethanolamines (NAPEs) and C20:1 MAG (MAG 20:1) reached the peak time in the sleep phase during DRF (Fig. 6e). While arrhythmic during DRF, myristoleoyl MAG (MAG14:1), oleoyl MAG (MAG 18:1) and N-acylethanolamines (NEA, endocannabinoids) as a class (represented by oleoyl ethanolamine/18:1-EA and palmitoyl ethanolamine/16:0-EA) exhibited robust daily oscillation in livers from NRF mice (Fig. 6e). This rhythmic pattern of precursors and components of endocannabinoids coexists with the rhythmicity of the fatty acid amide hydrolase FAAH (Fig. 3d).

Among the 489 mRNA-protein pairs, there is one gene (the gluconeogenic gene Pck1) that exhibited diurnal rhythmicity at four and the most levels of omics, followed by five genes (Plin5, Got1, Sun2, Hmgcr, and Abcg5) with diurnal activity across three levels of omics (Fig. 7b). Phase plot showed that almost all these genes exhibited various degrees of phase delay from a transcript-level rhythm to protein-level rhythms (Fig. 7b). This is matched with the patterns seen in circadian clock genes, such as Per2 and Nr1d1. Remarkably, the succinylated PCK1 rhythm exhibited an 8-h phase delay compared to its mRNA rhythm under NRF (Fig. 7b, c). The glutamic-oxaloacetic transaminase gene Got1 is an exception in that glycosylated GOT1 rhythm was 10.7 h advanced in phase compared to either protein or mRNA rhythms under DRF (Fig. 7b,c).

We analyzed the mRNA-protein pairs and found that the connections between diurnal multi-level proteomes and transcriptome exhibited a similar pattern between DRF and NRF (Fig. 7d). Transcript rhythms were mainly matched with rhythms in the unmodified protein or phosphorylated proteins, whereas unmodified protein rhythms were mainly matched with rhythms in mRNA, phosphorylated or ubiquitylated proteins. These findings revealed the complex network with regards to the biogenesis of diurnal rhythms in the liver.

Next, we compared the rhythmic pathways among different levels of omics, as measured by PSEA of rhythmic proteins and genes (Kuiper test, q < 0.05). Like with the gene-level comparison, most pathways had no matched rhythms across the omics under either DRF or NRF (Supplementary Fig. 7a). Connectivity maps for the rhythmic pathways (224 and 185 for DRF and NRF, respectively) revealed that strong connections under NRF including the mRNA-phosphorylation and phosphorylation-ubiquitylation links were replaced by connections such as the unmodified protein-glycosylation and phosphorylation-glycosylation links under DRF (Fig. 7e). Comparative analysis of rhythmic pathways between transcriptome and proteomes showed that non-coding RNA processing and metabolic process (ZT8-10) were found at the transcript and phosphorylated protein levels under NRF, while pathways related to translation (ZT18-20) were connected between transcript and protein levels under DRF (Supplementary Fig. 7b).

We next sought to validate the nutrient-responsive activity of PER2-pSer971. Ser971 is in the carboxyl-terminal region that is involved in the physical interaction with PPARγ (Fig. 7f). AML-12 mouse hepatocytes stably expressing enhanced green fluorescence protein (EGFP)/FLAG-tagged mouse PER2, as well as the non-phosphorylatable serine-to-alanine mutant (S971A), were generated by lentiviral particles. Using a polyclonal antibody raised against the pSer971 peptide of mouse PER2 (Supplementary Fig. 7c), we showed that PER2 Ser971 is endogenously phosphorylated in hepatocytes, the signals of which were abolished in PER2^S971A mutant (Fig. 7g). By immunoblotting analysis, we found that the levels of PER2-pSer971 increased as the concentration of extracellular glucose concentration raised from 0 to 5 mM (Fig. 7g). These data suggest that PER2-pSer971 is activated by nutrient flux in vivo.

Animal experiments were approved by the Laboratory Animal Welfare and Ethics Committee of Army Medical University (AMU), China (approval no. AMUWEC20201106). All experiments conformed to the relevant regulatory standards of AMU. Special pathogen-free (SPF) mice were purchased from Hunan SJT Laboratory Animal Co. and housed in a SPF barrier facility. C57BL/6J female mice at 7 weeks of age were grouped housed and entrained to a 12 h light:12 h dark cycle (light intensity 200 lux, humidity 40–60%, ambient temperature 21–23 °C) with free access to normal chow food (Jiangsu-Xietong, SFC9112-1010083, 22.8% kcal protein, 13.8% fat and 63.4% kcal carbohydrate, 3656 kcal/kg) and water.

C57BL/6J female mice at 7 weeks of age were grouped housed and entrained to a 12 h light:12 h dark cycle (light intensity 200 lux) with normal chow food and water ad libitum for at least 7 days before assigned to daytime-restricted feeding (DRF) and nighttime-restricted feeding (NRF) groups for 7 days. To synchronize the estrus cycle (4–5 days) of our animals, we co-housed these mice in a big rodent cage for 5–7 days before relocating to regular cages (n = 4 per cage) in our SPF facility. When assigning the mice to different groups, we randomized and controlled the body weight statistics to make sure the body weight is comparable between NRF and DRF groups. The DRF group had access to food for 12 h from ZT0 to ZT12 where ZT0 denotes light on at 9:00 am in China Standard Time (UTC + 8). The NRF group had access to food for 12 h from ZT12 to ZT24 (ZT0). At the end of time-restricted feeding, mice were euthanized by cervical dislocation and subjected to tissue collection every 4 h over 24 h (lipidomics) or 48 h (proteomics and PTM proteomics). Liver tissue was dissected (n = 4 per group per time-point totaling 56 or 96 samples). Tissues were snap frozen in liquid nitrogen and stored in a −80 °C fridge.

Mouse liver samples were pulverized in liquid nitrogen and sonicated in the urea lysis buffer containing 100 mM NH₄HCO₃ (pH 8), 8 M Urea and 0.2% Sodium dodecyl sulfate (SDS) for 5 min on ice. The lysate was centrifuged at 12,000 × g for 15 min at 4 °C, the supernatants were reduced by 10 mM DTT for 1 h at 56 °C, and subsequently alkylated with sufficient iodoacetamide (Sigma, #I6125-25G) for 1 h at room temperature in the dark. Next, samples were completely mixed with 4X volume of precooled acetone (Beijing Chemical Works, #11241203810051) by vortexing and incubated at −20 °C for at least 2 h. After centrifugation, the pellet was washed twice with cold acetone, and dissolved by dissolution buffer, which containing 0.1 M triethylammonium bicarbonate (TEAB, pH 8.5; Sigma, #T7408-500ML) and 6 M urea. Protein samples were subjected to 12% SDS-PAGE, all of which showed no visible sign of degradation.

120 μg liver proteins were aliquoted for peptide digestion by trypsin (Promega, #V5280), stopped with equal volume of 1% formic acid (HCOOH, Fisher Chemical, #A117-50), desalted in a C18 column and lyophilized. Lyophilized samples were dissolved in mobile phase A (2% acetonitrile or ACN, Fisher Chemical, #A955-4; adjusted to pH 10.0 by ammonium hydroxide) and fractionated in a C18 column (Waters BEH C18 4.6 × 250 mm, 5 μm) on a Rigol L3000 HPLC system at 50 °C (column oven) in an elution gradient mixed by mobile phase A and B (98% ACN, pH 10.0) at a flow rate of 1 mL/min as follows. 3–5% B, 10 min; 5–20% B, 20 min; 20–40% B, 18 min; 40–50% B, 2 min; 50–70% B, 3 min; 70–100% B, 1 min. Elutes were combined into 10 fractions, vacuum-dried and reconstituted in 0.1% formic acid in H₂O.

For data-dependent acquisition (DDA) spectrum library construction, shotgun proteomics analyses were performed using an EASY-nLC 1200 UHPLC system (Thermo Fisher, #LC140) coupled to an Q Exactive HF-X mass spectrometer (Thermo Fisher) via a nano-electrospray ion source (Nanospray Flex, spray voltage of 2.5 kV) at a capillary temperature of 320 °C. 1 μg sample containing iRT standard peptides (Biognosys, #Ki-3002-2) was analyzed by UHPLC at a flow rate of 600 nL/min in a linear gradient as follows. 5–8% B, 1 min; 8-30% B, 75 min; 30–50% B, 5 min; 50–95% B, 1 min; 95% B, 10 min; 95–5% B, 0.5 min; 5% B, 1 min; 5–95% B, 0.5 min; 95% B, 5 min; 95–5% B, 1 min. The Q Exactive HF-X was operated in Top40 mode with a full scan of 350–1500 m/z at a resolution of 120,000 (m/z 200). The automatic gain control (AGC) target value was set to 3 × 10⁶ at a maximum ion injection time of 80 ms. Precursor ions were fragmented by higher energy collisional dissociation (HCD) at a normalized collision energy (NCE) of 27% and analyzed in MS/MS, with the resolution set to 15,000 (m/z 200), an AGC of 5 × 10⁴, a maximum ion injection time of 45 ms, an intensity threshold of 1.1 × 10⁴ and a dynamic exclusion parameter of 20 s. For DIA acquisition, the m/z range was set to 350–1500 m/z at a resolution of 60,000 (m/z 200). The AGC was set to 5 × 10⁵ at a maximum ion injection time of 20 ms. Peptides were fragmented by HCD (NCE 27%) and fragment ions were analyzed in 50 scan windows in MS2, with the resolution set to 30,000 (200 m/z) and the AGC target value set to 1 × 10⁶.

DDA and DIA data were analyzed using Proteome Discoverer (Thermo, v2.4), Spectronaut (Biognosys, v9.0), and R against the UniProt mouse reference proteome (Proteome ID UP000000589, release 2019.01.18). DDA MS raw files were analyzed by Proteome Discoverer, and peak lists were searched against protein database. Cysteine carbamidomethylation (C/ + 57.021 Da) was set as a fixed modification; N-terminal acetylation (N-terminal/+42.011 Da) as N-Terminal modification; methionine oxidation (M/ + 15.995 Da) as dynamic modification. In peptide identification, maximum missed cleavage sites were set to 2, precursor mass tolerance set to 10 ppm and fragment mass tolerance set to 10 ppm. The false discovery rate (FDR) was set to 1% for both proteins and peptides, and was determined by searching a protein reverse database. MS1 label-free quantification (LFQ) was analyzed by maxLFQ algorithm and MS2 LFQ of DIA raw data was analyzed by Spectronaut with the following parameters, i.e., data extraction was set to “dynamic” (correction factor 1), identification set to “normal distribution p-value estimator” (q-value < 0.01), profiling strategy set to “iRT profiling” (q-value < 0.01), protein inference set to “from search engine”, protein quantity set to “average precursor quantity” and smallest quantitative unit set to “precursor ion” (summed fragment ions).

5 mg proteins were aliquoted for trypsin digestion, desalted in the C18 desalting column and lyophilized. Phospho-peptide fractions were enriched in an IMAC-Fe column (Fisher Scientific, #A32992) and incubated at room temperature for 30 min before centrifugation and lyophilization. Shotgun proteomics analyses were performed using an EASY-nLC 1200 UHPLC system (Thermo Fisher) coupled with an Q Exactive HF-X mass spectrometer (Thermo Fisher) operating in the DDA mode via a nano-electrospray ion source (Nanospray Flex, spray voltage of 2.3 kV) at a capillary temperature of 320 °C. 1 μg sample was analyzed by UHPLC at a flow rate of 600 nL/min in a 120 min linear gradient from 5 to 100% of eluent B (0.1% HCOOH in 80% ACN) in eluent A (0.1% HCOOH in H₂O) at a flow rate of 600 nL/min and detailed as follows, 5–10% B, 2 min; 10–30% B, 110 min; 30–50% B, 5 min; 50–95% B, 1 min; 95% B, 5 min. The Q Exactive HF-X was operated in Top30 mode with a full scan of 350–1500 m/z at a resolution of 120,000 (m/z 200). The AGC target value was set to 3 × 10⁶ at a maximum ion injection time of 80 ms. Precursor ions were fragmented by HCD at a NCE of 27% and analyzed in MS/MS, with the resolution set to 15,000 (m/z 200), an AGC of 5 × 10⁴, a maximum ion injection time of 100 ms, an intensity threshold of 5 × 10³ and a dynamic exclusion parameter of 30 s.

The resulting spectra were searched against the UniProt mouse reference proteome (Proteome ID UP000000589, release 2019.01.18) in Proteome Discoverer (v.2.2). Searched parameters: maximum missed cleavage sites of 2, mass tolerance of 10 ppm for precursor ion scans and a mass tolerance of 0.02 Da for the product ion scans. Cysteine carbamidomethylation (C/+57.021 Da) was set as a fixed modification; N-terminal acetylation (N-terminal/+42.011 Da) as N-Terminal modification; phosphorylation of serine, threonine and tyrosine residue (S, T, Y/+79.966 Da) and methionine oxidation (M/+15.995 Da) as dynamic modifications; FDR < 1% at levels of PSM and protein (matched with at least one unique peptide).

10 mg proteins were aliquoted for trypsin digestion, and incubated with anti-Ubiquitin Remnant Motif (K-ε-GG) beads (CST #5562) in MOPS IAP buffer (50 mM MOPS, 10 mM KH₂PO₄, and 50 mM NaCl) for 2 h at 4 °C. Shotgun proteomics analyses were performed using an EASY-nLC 1200 UHPLC system (Thermo Fisher) coupled with an Q Exactive HF-X mass spectrometer (Thermo Fisher) operating in the DDA mode via a nano-electrospray ion source (Nanospray Flex, spray voltage of 2.3 kV) at a capillary temperature of 320 °C. 1 μg sample was analyzed by UHPLC at a flow rate of 600 nL/min in a 120 min linear gradient from 5 to 100% of eluent B (0.1% HCOOH in 80% ACN) in eluent A (0.1% HCOOH in H₂O) at a flow rate of 600 nL/min and detailed as follows, 5–10% B, 2 min; 10–30% B, 110 min; 30–50% B, 5 min; 50–95% B, 1 min; 95% B, 5 min. The Q Exactive HF-X was operated in Top30 mode with a full scan of 350–1500 m/z at a resolution of 120,000 (200 m/z) with an AGC target value of 3 × 10⁶ and a maximum ion injection time of 80 ms. Precursor ions were fragmented by HCD at a NCE of 27% and analyzed in MS/MS, with the resolution set to 15,000 (m/z 200), an AGC of 5 × 10⁴, a maximum ion injection time of 100 ms, an intensity threshold of 5 × 10³ and a dynamic exclusion parameter of 30 s.

The resulting spectra were searched against the UniProt mouse reference proteome (Proteome ID UP000000589, release 2020.01.08) in Proteome Discoverer (v.2.4). Searched parameters: maximum missed cleavage sites of 2, mass tolerance of 10 ppm for precursor ion scans and a mass tolerance of 0.02 Da for the product ion scans. Cysteine carbamidomethylation (C/ + 57.021 Da) was set as a fixed modification; N-terminal acetylation (N-terminal/+42.011 Da) as N-Terminal modification; GG of lysine (K/ + 114.043 Da) and methionine oxidation (M/ + 15.995 Da) as dynamic modification; FDR < 1% at levels of PSM and protein (matched with at least one unique peptide).

10 mg proteins were aliquoted for trypsin digestion, and incubated with anti-succinyl-lysine motif beads (CST #13764) in MOPS IAP buffer for 2 h at 4 °C. Shotgun proteomics analyses were performed using an EASY-nLC 1200 UHPLC system (Thermo Fisher) coupled with an Q Exactive HF-X mass spectrometer (Thermo Fisher) operating in the DDA mode via a nano-electrospray ion source (Nanospray Flex, spray voltage of 2.3 kV) at a capillary temperature of 320 °C. 1 μg sample was analyzed by UHPLC at a flow rate of 600 nL/min in a 120 min linear gradient from 5 to 100% of eluent B (0.1% HCOOH in 80% ACN) in eluent A (0.1% HCOOH in H₂O) at a flow rate of 600 nL/min and detailed as follows, 5–10% B, 2 min; 10–30% B, 110 min; 30–50% B, 5 min; 50–95% B, 1 min; 95% B, 5 min. The Q Exactive HF-X was operated in Top30 mode with a full scan of 350–1500 m/z were acquired at a resolution of 120,000 (200 m/z) with an AGC target value of 3 × 10⁶ and a maximum ion injection time of 80 ms. Precursor ions were fragmented by HCD at a NCE of 27% and analyzed in MS/MS, with the resolution set to 15,000 (200 m/z) with an AGC target value of 5 × 10⁴, a maximum ion injection time of 100 ms, an intensity threshold of 5 × 10³ and the dynamic exclusion parameter of 30 s.

The resulting spectra were searched against the UniProt mouse reference proteome (Proteome ID UP000000589, release 2020.01.08) in Proteome Discoverer (v.2.4). Searched parameters: maximum missed cleavage sites of 2, mass tolerance of 10 ppm for precursor ion scans and a mass tolerance of 0.02 Da for the product ion scans. Cysteine carbamidomethylation (C/ + 57.021 Da) was set as a fixed modification; N-terminal acetylation (N-terminal/+42.011 Da) as N-Terminal modification; succinylation of lysine (K/ + 100.016 Da) and methionine oxidation (M/ + 15.995 Da) as dynamic modification; FDR < 1% at levels of PSM and protein (matched with at least one unique peptide).

500 μg proteins were aliquoted for trypsin digestion, desalted, and vacuum-dried. Peptides were dissolved in 80% ACN (0.1% trifluoracetic acid) and enriched for N-glycosylated peptides in the HILIC column (Bonna-Agela, #VH950010). The mix was first washed in 0.1% trifluoracetic acid, then in 50 mM NH₄HCO₃-H₂¹⁸O and digested in PNGase-F (Sigma, #P7367) overnight at room temperature. Shotgun proteomics analyses were performed using an EASY-nLC 1200 UHPLC system (Thermo Fisher) coupled with an Q Exactive HF-X mass spectrometer (Thermo Fisher) operating in the DDA mode via a nano-electrospray ion source (Nanospray Flex, spray voltage of 2.3 kV) at a capillary temperature of 320 °C. 1 μg sample was analyzed by UHPLC at a flow rate of 600 nL/min in a 120 min linear gradient from 5 to 100% of eluent B (0.1% HCOOH in 80% ACN) in eluent A (0.1% HCOOH in H₂O) at a flow rate of 600 nL/min and detailed as follows, 5–10% B, 2 min; 10–40% B, 105 min; 40–50% B, 5 min; 50–90% B, 3 min; 90–100% B, 5 min. The Q Exactive HF-X was operated in Top40 mode with a full scan of 350–1500 m/z with a resolution of 60,000 (200 m/z) with an AGC target value of 3 × 10⁶ and a maximum ion injection time of 20 ms. Precursor ions were fragmented by HCD at a NCE of 27% and analyzed in MS/MS, with the resolution set to 15,000 (200 m/z) with an AGC target value of 5 × 10⁴, a maximum ion injection time of 80 ms, an intensity threshold of 1.3 × 10⁴, and the dynamic exclusion parameter of 30 s.

The resulting spectra were searched against the UniProt mouse reference proteome (Proteome ID UP000000589, release 2019.01.18) in Proteome Discoverer (v.2.2). Searched parameters: maximum missed cleavage sites of 2, mass tolerance of 10 ppm for precursor ion scans and a mass tolerance of 0.02 Da for the product ion scans. Cysteine carbamidomethylation (C/ + 57.021 Da) was set as a fixed modification; N-terminal acetylation (N-terminal/+42.011 Da) as N-Terminal modification; glycosylation of asparagine (deamidation ¹⁸O (N)/ + 2.988 Da) and methionine oxidation (M/ + 15.995 Da) as dynamic modification; FDR < 1% at levels of PSM and protein.

Lipids were extracted from approximately 20 mg tissues using a modified version of the Bligh and Dyer’s method³⁴. Liver tissues were homogenized in 750 µL of chloroform:methanol 1:2 (v/v) with 10% deionized water on a bead ruptor (OMNI, USA). The homogenate was then incubated at 1500 rpm for 1 h at 4 °C. At the end of the incubation, 350 µL of deionized water and 250 µL of chloroform were added to induce phase separation. The samples were then centrifuged and the lower organic phase containing lipids was extracted into a clean tube. Lipid extraction was repeated once by adding 500 µL of chloroform to the remaining tissues in aqueous phase, and the lipid extracts were pooled into a single tube and dried in the SpeedVac under OH mode. Polar lipids were analyzed using an Exion-UPLC system coupled with a triple quadrupole/ion trap mass spectrometer (6500 Plus Qtrap; SCIEX)^35,66. Separation of individual lipid classes of polar lipids by normal phase (NP)-HPLC was carried out using a Phenomenex Luna 3 µm-silica column (internal diameter 150 × 2.0 mm) with the following conditions: mobile phase A (chloroform:methanol:ammonium hydroxide, 89.5:10:0.5) and mobile phase B (chloroform:methanol:ammonium hydroxide:water, 55:39:0.5:5.5). MRM transitions were set up for comparative analysis of various polar lipids. Glycerol lipids including diacylglycerols and triacylglycerols were quantified using a modified version of reverse phase HPLC/MRM. Separation of neutral lipids were achieved on a Phenomenex Kinetex-C18 2.6 µm column (4.6 × 100 mm) using an isocratic mobile phase containing chloroform:methanol:0.1 M ammonium acetate 100:100:4 (v/v/v) at a flow rate of 170 µL for 17 min.

N-acylethanolamines (NAEs) and N-acyl phosphatidylethanolamines (NAPEs) were extracted from liver tissues using a modified method of Bligh and Dyer’s extraction⁶⁷. The organic phase was extracted and dried in SpeedVac under organic mode. Samples were resuspended in acetonitrile: isopropanol (1:2 v/v) containing appropriate concentrations of internal standards including d7-PE33:1 from Avanti Polar Lipids, d8-20:4-EA, d4-16:0-EA, d4-22:6-EA and d5-MAG-20:4 from Cayman Chemicals. Samples were analyzed on a ThermoFisher U3000 DGLC coupled to Sciex 6500 Plus QTRAP under the electrospray ionization mode. Individual lipids were separated on an Agilent Zorbax Eclipse Plus column (100 × 2.1 mm, 1.8 µm) using Mobile Phase A (10 mM ammonium formate: acetonitrile: isopropanol 50:30:20, pH 8) and Mobile Phase B (10 mM ammonium formate: acetonitrile: isopropanol 1:9:90, pH 8). NAEs and NAPEs were quantitated by referencing to spiked internal standards.

Integrative network analysis of the hepatic phosphoproteome, proteome, and ubiquitylome was performed using the DIABLO (Data Integration Analysis for Biomarker discovery using Latent variable approaches for Omics studies) model in the R package mixOmics⁴⁷. First, omics features were selected using sparse Partial Least Squares Discriminant Analysis (sPLS-DA) based on three subsets of input data, i.e., morning (ZT0-4), evening (ZT8-12) or night (ZT16-20). Abundance profiles of 2139 rhythmic phosphoproteins, 1983 rhythmic proteins, and 483 rhythmic ubiquitylated proteins found in livers from either DRF or NRF mice served as input data. The evening subset was chosen for modeling study for its low BER. Specifically, the following modeling parameters associated with the evening subset were used after tuning: ncomp = 3 (number of model components); the number of rhythmic proteins to consider: 20; the number of rhythmic phospho-proteins to consider: 6; the number of rhythmic ubiquitylated proteins to consider: 5, 5. Data were then modeled in the DIABLO framework. Plots of experimental groups (function: plotIndiv) provide a visual representation of samples in the DIABLO framework. Correlation Circle plot (function: circosPlot), Clustered Image Maps for DIABLO (function: cimDiable), and plots of variable loadings (importance of the feature/molecule) were produced in the mixOmics R package.

Alpha mouse liver 12 (AML-12) cell line was established from mouse hepatocytes (source: a 3-month-old CD1 strain MT42 mouse line) that ectopically express human transforming growth factor alpha (ATCC #CRL-2254). AML12 cells were authenticated by their general morphology, short tandem repeat microsatellite marker analysis, capacity to form lipid droplets upon oleic acid treatment and RT-qPCR for mouse genes. AML-12 cells were cultured at 37 °C and 5% CO₂ in DMEM/F12 (Procell #PM150312) supplemented with 10% fetal bovine serum (Procell #164210-50), 10 µg/mL insulin, 5.5 µg/mL transferrin, 5 ng/mL selenium, 40 ng/mL dexamethasone and 1% penicillin/streptomycin (100 U/mL and 0.1 mg/mL, respectively; Procell #PB180120↗), which is formulated as the AML-12 culture medium (Procell #CM-0602). Cells were seeded to 6-well plates (NEST #703001) or 10-cm dish (NEST #704002). When reaching 80-90% confluence, cells were treated in DMEM (Gibco #11966025) containing various concentrations of glucose (BioReagent grade, Sigma #G7021) for 8 h.

Mouse Per2 cDNA (NM_011066.3↗, CDS 175..3948 nucleotides) or its S971A variant was generated by a Gibson Assembly Cloning Kit (NEB #E5510) and cloned into the lentiviral vector pLenti-hEF1a-PuroR-CMV-MCS-EGFP-3xFlag (BamHI/XmaI). Lentiviral particles, i.e., Lenti-hEF1a-PuroR-CMV-Per2-EGFP-3xFlag (3.44E + 08 transduction unit per mL or TU/mL) and Lenti-hEF1a-PuroR-CMV-Per2^S971A-EGFP-3xFlag (3.76E + 08 TU/mL) were prepared by the vendor (Shanghai Tailtool Bioscience Co. Ltd). AML-12 mouse hepatocytes in early passages were seeded at 10,000 cells per well into a 6-well cell culture microplate and transduced with 5 µg/mL polybrene and lentiviral particles at a multiplicity of infection (MOI) of 50 over the night. These cells express PER2 or PER2^S971A as a fusion protein that is tagged in the carboxyl-terminal with EGFP and 3xFLAG peptide. Stable cells were cultured at 5 µg/mL puromycin (Beyotime #ST551) for 7–10 days and selected for EGFP biofluorescence under a Nikon-U fluorescence microscope.

Cells were lysed in RIPA buffer (ThermoFisher #89901) supplemented with proteinase inhibitor cocktail (Roche #5892953001) and protein phosphatase inhibitor cocktail (PhosSTOP, Roche #4906845001). 15 µg of protein lysates were electrophoresed on SDS-PAGE gels at 110 V and transferred to activated PVDF membranes (Bio-Rad #1620177). Block the membranes in 5% skim milk or 5% BSA solution (PER2-pSer971) for 1.5 h and incubated with primary antibodies at a dilution of 1:1000 (except β-Actin 1:50,000) and 4 °C overnight. PER2 rabbit polyclonal antibody (pAb) (Abclonal #A13168, RRID:AB_2760019, 1:1000); PER2-pSer971 rabbit pAb (Abclonal #WG-00762P, 1:1000); β-Actin rabbit mAb (Abclonal #AC026, RRID:AB_2768234, 1:50,000); DYKDDDDK (FLAG) tag rabbit pAb (Proteintech, 20543-1-AP, RRID: AB_11232216, 1:1000). Immunoblots were visualized by peroxidase conjugated secondary antibodies, i.e., HRP-linked anti-mouse IgG (Cell Signaling Technology, Cat# 7076, RRID: AB_330924, 1:5000) and HRP-linked anti-rabbit IgG (Cell Signaling Technology, Cat# 7074, RRID: AB_2099233, 1:5000), and enhanced chemiluminescent substrate (Bio-Rad #1705062) in the Azure C500 imaging system (see Fig. 7g and associated Source Data file). PVDF membranes were stripped in antibody removal solution (Beyotime #P0025B) for 20 min and washed three times in TBST before re-starting the blocking and immunoblotting procedure. Experiments were reproduced with similar results (n = 8 biological replicates from 4 independent experiments). PER2-pSer971 rabbity polyclonal antibody was generated by Abclonal Inc. Phosphopeptide GRA(S-p)PPLFQSR-C (mPER2 968-978 AA) was conjugated to KLH and immunized to SPF rabbits (Oryctolagus cuniculus). The antibody was purified from sera by affinity purification and screened for positive signals against the phosphopeptides and negative signals against the control peptide GRASPPLFQSR-C.

Phase set enrichment analysis was performed on oscillating genes in different tissues to identify phase-clustered pathways via PSEA1.1_VectorGraphics.jar⁶⁸. Source files (v7.1) are obtained from the Molecular Signatures Database (Broad Institute, MA, USA). Mouse gene nomenclature was converted to human gene nomenclature via “Human and Mouse Homology Classes with Sequence information” (curated by http://www.informatics.jax.org↗). Parameters: domain 0 to 24, minimal number 10, and maximal number 10000. Q-value < 0.05 is considered as statistical significance. Results were sorted by Vector-average magnitude and Vector-average value, respectively before visualization.

Sample size was determined by Guidelines for Genome-Scale Analysis of Biological Rhythms⁶⁹. Rhythmic proteins, sites or lipids were considered as reaching statistical significance in three rhythmicity methods, i.e., meta2d_pvalue < 0.05 (MetaCycle)⁴⁶, adjusted P-value < 0.05 (RAIN)⁴⁴ and “Presence of rhythmicity p-value < 0.05” (Circacompare)⁴⁵. Two-sided unpaired Student’s t-test with or without Bonferroni correction and two-sided unpaired Wilcoxon test were performed in R (version 4.1.3)/RStudio (2023.03.0 Build 386) as indicated, and P-value < 0.05 is considered as statistical significance. Heatmap of diurnal expression profile was visualized by R package: pheatmap (version 1.0.12). Data were plotted via R packages: ggplot2 (version 3.3.6) and ggpubr (version 0.4.0). Venn and chord diagrams were plotted via R package: VennDiagram (version 1.7.3) and Circlize (version 0.4.15).

Further information on research design is available in thelinked to this article. Nature Portfolio Reporting Summary

Supplementary Information Reporting Summary Peer Review File Supplementary Table 1

Source Data