In-Vivo Quantitative Proteomics Reveals a Key Contribution of Post-Transcriptional Mechanisms to the Circadian Regulation of Liver Metabolism

Liver samples were harvested from wildtype C57BL/6 mice kept one day in constant darkness after being entrained to a 12–12 h light-dark schedule. Protein lysates were prepared from the mouse livers, which were collected at intervals of 3 h over two circadian cycles (4 mice per time point; n = 64 total mice). For each time point, protein extracts from the four livers were mixed in equal amounts to have a single sample per time point. We decided to use SILAC in an in vivo format for the proteomic quantification method [23], [24]. An internal, spike-in standard mix was constructed by combining equal amounts of protein lysates from liver samples of two heavy SILAC labeled mice collected in anti-phase (see Material and Methods). The pooled lysates of each of the 16 time points was mixed with equal amounts of the internal standard prior to digestion. Resulting peptides were then separated into six fractions and measured on a linear ion trap – Orbitrap mass spectrometer (Figure 1A). The experiment was performed in technical triplicates resulting in 288 liquid chromatography (LC)-MS/MS files that were subsequently processed in the MaxQuant software environment [25]. The relative abundance of the liver proteome is calculated for each time point by taking the ratio of the mass spectrometric signal for individual proteins to the signal of the spiked in heavy SILAC standard. For circadian analysis we filtered out those protein groups with accurate quantification values in less than half of the samples measured; the resulting dataset contained 3132 protein groups (Table S1).

Reproducibility was tested by comparing the measurements to each other and calculated Pearson correlation coefficients. For all of the comparisons we obtained high r values (between 0.6 and 0.93) as illustrated for technical and biological replicates in Figure 1B. The stringently quantified proteome dataset of 3132 contains proteins from a broad range of metabolic and cellular processes. Coverage of transcription factors and cell cycle proteins in these post-mitotic cells was relatively low, as expected. We did not obtain quantification values for clock proteins, likely due to their low abundance particularly at some times of the cycle. Therefore proper entrainment of the mice was confirmed by assessing the expression profile of Bmal1 and Per2 mRNA as well as of PER2 protein in the collected liver samples (Figure S1). Nevertheless, the quantified liver proteome covered extensively many liver specific pathways like fatty acid and drug metabolism (Figure 1C). In addition, more than 50% of coverage was observed for more general cellular processes and components such the ribosome, proteasome and spliceosome (Figure 1C).

To identify proteins with circadian rhythms of abundance in the mouse liver quantified dataset we adapted a statistical algorithm (Materials and Methods). Specifically, we determined goodness-of-fit of the expression ratios to a cosine curve with a period of 23.6 hours, the circadian period reported for the used mouse strain [26]. To calculate the rate of false discovery the experimental data were repeatedly scrambled and fitted to the cosine curve, preserving the technical triplicates in the randomization. This allows an estimate of the frequency that a false observation matches the curve by chance. With this statistical analysis, at a false discovery rate (FDR or q-value) <0.33, we identified 186 proteins (Table S2) from the total 3132 dataset (Table S1) that showed circadian rhythmic profiles of abundance in the mouse liver across two consecutive cycles. Based on this analysis, we thus estimate that at least 6% of our mouse liver proteome dataset shows circadian oscillations. We additionally compared the statistical analysis performed by our method to the JTK Cycle method, a standard method used in the circadian field to detect rhythmicity [27]. The results indicated an excellent correlation (Pearson r>0.9) between the two methods for q-value estimation in the total quantified dataset (Figure S1C) as well as phase determination of statistical significant cycling proteins (Figure S1D). This shows that the statistical algorithm to detect oscillations included in the Perseus software package may be widely applicable in circadian studies.

We assessed the total abundance of each protein in the samples by using the added peptide intensity obtained in the MS analysis normalized to their molecular weight. The abundance profile of the cycling proteome was very similar to the one observed for the entire quantified proteome (Figure 2A). This indicates that cycling proteins identified in our analysis are not biased against low abundance. Although the amplitude of oscillations has been studied extensively at the transcript level, there is currently no information about the accurate fold change of the cycling proteome. Using the logarithmic normalized expression ratios obtained for every data point, we calculated the amplitude of the abundance for each oscillating protein across the circadian day. The mean of the distribution of the total fold change for rhythmic proteins in the mouse liver is 1.38; therefore the large majority of rhythmic proteins change less than two fold (Figure 2B). The abundance of cycling proteins does not correlate with the coefficient of variation of the oscillation (Pearson r = 0.018), indicating that circadian rhythms were equally detectable in high and low abundance proteins. Next we tested whether there were categorical protein annotations significantly different from the overall fold change distribution by applying a recently developed annotation enrichment algorithm [28]. Proteins cycling with small amplitude are statistically enriched in those annotated as being modified by acetylation as well as ubiquitin-like modifier proteins whereas glycosylated proteins show cycles with large amplitudes (Benjamini Hochberg FDR<0.05 for all; Figure S1E).

Enrichment analysis of the circadian proteome compared to the total quantified proteome revealed that several essential liver metabolic Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways such as drug and bile acids metabolism were over-represented (Fisher test p<0.05) (Figure 2C, Table S3). Interestingly, a significant fraction of the membrane bound proteome (Figure 2D, Table S3) is under circadian regulation in the mouse liver. In addition, secreted proteins from both extracellular matrix as well as plasma- as those involved in blood coagulation- all synthesized in the liver, are also enriched (Fisher test p<0.01) in the circadian proteome suggesting that the circadian clock in the liver may play a role in the generation of rhythms of plasma and extracellular proteins (Figure 2C and 2D, Table S3).

Taken together, we report here the most comprehensive and accurately measured circadian proteome of the mouse liver to date and find that it is significantly enriched in essential protein categories whose temporal regulation has not been previously documented.

We performed phase dependent hierarchical cluster analysis of the mouse liver circadian proteome using the normalized logarithmic expression ratios for all sampled circadian times (CT). The heat map representation of the ratios showed cycling proteins with phases distributed across the two consecutive days (Figure 3A). The overall pattern of these phases was similar to previously reported circadian profiles of transcripts in this organ [3]–[6]. The analysis resulted in two major branches in the dendrogram which mostly segregate proteins peaking during the day or during the night (Figure 3A, 3B). Almost 2/3 of the cycling proteome displayed night phases (CT12 to CT24) and 1/3 day phases (CT0 to CT12). Mice are nocturnal animals hence behavioral and metabolic activity is predominant during the dark phase, which may explain the larger number of night peaking proteins in the liver. We next asked if any protein category was statistically different in these two main clusters, using a Fisher exact test cut-off of p<0.02. Secreted and extracellular proteins tended to be rhythmic with abundance peak during the day while proteins associated to membrane, to the endoplasmic reticulum (ER) as well as to the Golgi apparatus mainly peaked at night. Moreover, proteins from complement and coagulation cascades oscillated with exclusively day phases (Figure 3B and S2A) while many rhythmic proteins with nocturnal phases are involved in drug metabolism, bile secretion and protein processing in ER (Figure 3B and S2B).

To identify additional phase dependent enriched categories of cycling proteins in the mouse liver we directly tested for phase enrichment of categorical annotations in a wider set (838) of cycling candidates using a loser q-value cut-off (<0.66). To determine the significance of the cycling annotation distribution test, we employed a Fisher Exact Test cut-off of FDR<0.02, resulting only in categories of proteins with phases highly significant clustered at specific time of the circadian cycle. Day-enriched Gene Ontology Cellular Component (GOCC) protein categories comprised mainly nuclear and extracellular proteins while mitochondrial associated proteins peaked at the day-night interphase (Figure S2C). In contrast, proteins associated to Golgi displayed statistically enriched night phases. The cycling annotation distribution analysis also highlighted KEGG metabolic pathways enriched at specific times of the circadian cycle. Drug metabolism and protein processing in ER were remarkably enriched with nocturnal phases which agree with the exclusive presence of proteins from these pathways in the night hierarchical cluster described above (Figure 3A, S2B and S2D). The complement and coagulation cascades showed, in contrast, day-enriched phases as mentioned above too (Figure 3A, S2A and S2D).

Together our enrichment analyses revealed that a broad range of cellular and metabolic components are subjected to temporal regulation by means of protein abundance. A considerable number of cycling proteins belong to categories that have not previously been described to be under circadian control in the mouse liver.

Circadian rhythms of mRNA in the mouse liver have been widely studied, however it still remains unresolved whether and how the temporal changes in transcripts translate to global oscillations of protein abundance. To this end, we compared our circadian proteome to a microarray study of the mouse liver transcriptome from the Hogenesch group [3]. That experiment was performed with 1 h resolution and around 10% of the transcriptome was found to be cycling, depending on the statistical algorithm employed. We matched each protein group of our quantified dataset to its corresponding Affymetrix entries which gave us a final dataset of 3046 protein groups. Using the same statistical algorithm that we applied to the proteome we assessed both circadian oscillations for proteins and transcripts. By using a cut-off of q-value<0.33 we identified 181 protein groups (6%) with significant circadian rhythmicity, among them 151 also showed rhythms at the mRNA level (Figure 4A, Table S4). Comparison of q-value distributions for the total dataset of protein and mRNA indicates that most of the cycling proteins with arrhythmic mRNA showed q-values for their transcripts far from the applied cut-off (Figure S3A). Thus, around 20% of cycling proteins do not seem to be accompanied by statistical significant oscillations of their mRNA based on the gene array data. Circadian gene expression can also be assessed by analyzing cycling binding patterns of CLOCK and BMAL transcription factors in promoters. Therefore, in addition to this gene expression study we used a recently reported set of BMAL1 target genes found by chromatin immunoprecipitation (ChIP) [29]. In that data set, we found only four additional genes where the corresponding protein was cycling but the transcript did not show statistically significant oscillation. Together our statistical analyses suggest that around 20% of the cycling proteins do not show circadian regulation of their corresponding transcript, implying wide-spread circadian post-transcriptional control of protein abundance.

Given that very recent evidence indicated that post-transcriptional mechanisms appear to largely determine the phase of mRNA oscillations [16], [30]–[32], we globally explored the distribution of phases for oscillating transcripts and proteins. We found that phases of rhythmic transcripts matching cycling proteins in the mouse liver were evenly distributed throughout the cycle with higher frequency during the night (Figure 4B upper graph) similarly to what was previously reported for the total circadian liver transcriptome [3]–[6]. In contrast, the phases of cycling proteins with rhythmic transcripts were distributed in two main clusters, a smaller one centered at the middle of the day and a larger one in the second part of the night. Surprisingly, very few proteins peaked in the first hours of the day (Figure 4B lower graph). The divergence in the distribution of transcript and protein phase suggests that the cycles of protein abundance in the mouse liver do not necessarily reflect mRNA changes and instead are influenced by post-transcriptional mechanisms. For the clock genes themselves, a characteristic time delay (usually 4–6 h) between mRNA and protein expression has been described [17]. However, it is poorly understood whether this is a general feature of all cycling genes and proteins. To address this question using the quantitative nature of our circadian proteome data, we calculated the time delay between the peak of expression of each rhythmic transcript and its oscillating protein. The total distribution of time delays showed with a preferential window of 2 to 6 h between transcript and protein peaks for almost 50% of the cycling mRNA and proteins (Figure S3C). Strikingly, around 40% of oscillating proteins peaked more than 6 h later than their corresponding transcripts. This data indicates that regulation of time delay between peaks of mRNA and protein is an important feature of circadian biology, implying general post-transcriptional mechanisms that define the overall phase-tuning of the proteome.

Principal component analysis (PCA) of the cycling proteome and transcriptome showed excellent clustering of the technical triplicates and biological replicates in the two-dimensional graph (Figure 4C). The graphical representation of the data on the basis of the two main PCA components resembles an analog clock (Figure 4D), indicating that time is the main component accounting for the overall difference between both datasets. We then calculated the angles corresponding to the median value for each time point experimentally analyzed in the proteome and for the published transcriptome. These angles directly visualize (Figure 4D) the characteristic time delay between mRNA and protein cycles mentioned above (Figure S3C). Very interestingly, the length of this time lag varies across the circadian cycle (Figure 4D). Together the comparison of the phases of the transcriptome vs. the proteome clearly demonstrates that circadian protein oscillations are shaped post-transcriptionally.

Having established that the circadian clock drives liver metabolism not only at the level of mRNA but also by precise regulation of the phase of protein abundance, we next examined individual metabolic processes. We found that rhythmic proteins associated to specific functions oscillated with similar phases regardless of whether their transcripts were cycling or, if so, when they were expressed. This implies that circadian post-transcriptional regulation coordinates individual metabolic pathways – something that became especially obvious for the metabolism of xenobiotics. Phases of abundance of crucial components of this pathway tightly clustered at the end of the night while there was no obvious coordination in the phases of the cycling mRNAs (Figure 5A). Moreover, temporal coordination of protein abundance of drug metabolism in the liver was not restricted to enzymes involved in the three phases of the detoxification mechanism. Strikingly, it also extended to membrane transporters responsible for the intake of xenobiotic substances from the blood circulation into the liver and/or their subsequent secretion into bile after being conjugated (Figure 5A and 5B). Supporting a metabolic function of this coordinated expression, we found that the abundance profiles of rhythmic enzymes involved in detoxification strongly correlated with the oscillations of xenobiotics levels reported recently in a liver circadian metabolome study [10] (Figure 5C).

A substantial number of proteins involved in complement and coagulation pathways oscillate in abundance, with their phases significantly clustered during the day. Their phases of abundance did not uniformly reflect the cycling patterns of their corresponding transcripts (). Figure S3D

Temporal regulation of gene expression mediated by the circadian transcription factor heat shock factor 1 (HSF1) has been reported for some heat shock proteins [33], however, our study shows for the first time that the protein folding machinery in the mouse liver is also regulated in a circadian manner at the protein level. We found many chaperones undergoing daily oscillations of abundance in the mouse liver with phases remarkably clustered at night. Heat shock proteins present an interesting contrast to the examples described above, in that protein cycles of most of them closely mirror their transcriptional profiles (Figure S3E) showing short time delays between mRNA and protein peaks. Thus, circadian regulation of protein abundance for this pathway seems to be extensively determined at the level of mRNA in contrast to what we observed for the metabolism of xenobiotics and coagulation cascades.

In summary, our analyses indicate that post-transcriptional mechanisms have a large but diverse influence on the oscillation of the proteome essential for different cellular processes. Specifically, the contributions of this circadian post-transcriptional control appear to differ among diverse metabolic pathways.

All mice were bred and maintained in the animal facility of the Max Planck Institute of Biochemistry according to institutional guidelines and all animal experiments were approved by the government agencies of Oberbayern. Eight-week-old C57BL/6 mice were house in light-tight boxes with free access to food and water and entrained to a 12–12 h light-dark schedule for ten days before being transfer to complete darkness. After one day in constant darkness, mice were sacrificed at 3 h intervals over two days. Prior to liver excision, mice were perfused with ice-cold PBS to remove blood content. Livers were then quickly frozen in liquid nitrogen followed by storage at −80°C. SILAC mice [23] were kept in the same conditions; two animals were sacrificed in anti-phase at CT3 and CT15, respectively, and livers excised as described above.

Protein extracts from mouse livers were obtained as previously described [56], [57]. Briefly, homogenization of 1 mg of liver was done in 1 ml of 0.1 mM Tris-HCl pH 7.6 supplemented with complete protease and phosphatase inhibitor cocktails (Roche) using an Ultra Turbax blender (IKA) at maximum speed at 4°C for 30–60 seconds. Sodium dodecyl sulfate (SDS) and dithiothreitol (DTT) were added to the homogenates to a final concentration of 4% and 0.1 mM, respectively, followed by brief sonication to reduce viscosity. After 5 min incubation at 95°C the mixture was then cleared by centrifugation at 16,000× g at room temperature for 10 min. Protein content was determined by measurements of tryptophan fluorescence as previously described [56]. Sixteen protein extract pools, corresponding to the samples collected at different circadian times, were obtained by mixing equal amounts of protein liver extracts from each of the four mice sacrificed at any given time point. Similarly the SILAC protein liver mix was obtained by adding equal amounts of the protein extracts obtained from the two SILAC liver samples collected in anti-phase.

For sample preparation we used 100 µg of protein extract from each circadian time pool mixed with 100 µg of SILAC protein pool. The protein mixes were concentrated in 30 k Microcon filtration devices (Millipore) to a final volume of 30 µl and then processed by the FASP procedure [57]. Briefly, the samples were mixed with 0.2 ml of 8 M urea in 0.1 M Tris/HCl pH 8.5 (UA), loaded into 30 k Microcon filtration devices (Millipore) and centrifuged at 14,000× g for 15 min. The concentrates were diluted in the devices with 0.2 ml of UA solution and centrifuged again. After centrifugation the concentrates were mixed with 0.1 ml of 50 mM iodoacetamide in UA solution and incubated in the dark at room temperature for 30 min. After centrifugation for 15 min the concentrate was diluted with 0.2 ml UA solution and concentrated again by centrifugation. This step was repeated twice. Next, the concentrate was diluted with 0.1 ml of 40 mM NaHCO₃ and concentrated again twice. Subsequently, 2 µg of Lysyl Endopeptidase (Wako Chemicals) in 40 µl of 40 mM NaHCO3 was added to the filter and the samples were incubated at room temperature overnight. The peptides were collected by centrifugation of the filter followed by two additional 30 µl washes with 40 mM NaHCO₃. The concentration of peptides was determined by UV-spectrometry using an extinction coefficient of 1.1 for 0.1% (g/l) solution at 280 nm.

Peptides were dissolved in 200 µL in Britton & Robinson buffer composed of 20 mM CH₃COOH, 20 mM H₃PO₄, and 20 mM H₃BO₃, and NaOH, pH 11. The peptides were separated by a pipette-based anion exchanger method [57]. Briefly, the pipette based column was assemble by stacking 6 layers of a 3M Empore Anion Exchange disk (Varian, 1214-5012) into a 200 µl micropipette tip. For column equilibration and elution of fractions Britton & Robinson buffer titrated with NaOH to the desired pH was used. Peptides were loaded at pH 11 and fractions were subsequently eluted with buffer solutions of pH 8, 6, 5, 4, and 3, respectively.

All mass spectrometric (MS) experiments were performed on a nanoflow HPLC system (Proxeon Biosystems, now Thermo Fisher Scientific) connected to a hybrid LTQ-Orbitrap (Thermo Fisher Scientific, Bremen, Germany), equipped with a nanoelectrospray ion source (Proxeon Biosystems, now Thermo Fisher Scientific). Peptide mixtures were separated by reversed phase chromatography using in-house-made C₁₈ microcolumns with a diameter of 75 µm packed with ReproSil-Pur C18-AQ 3-µm resin (Dr. Maisch GmbH, Ammerbuch-Entringen, Germany) in 4 hours LC gradient from 3% to 75% acetonitrile in 0.5% acetic acid at a flow rate of 200 nl/min and directly electrosprayed into the mass spectrometer. The LTQ-Orbitrap was operated in the positive mode to simultaneously measure full scan MS spectra (from m/z 300–1650) in the Orbitrap analyzer at resolution R = 60 000 following isolation and fragmentation of the ten most intense ions in the LTQ part by collision-induced dissociation. Raw MS files were processed with MaxQuant (version. 1.1.1.9), a freely available software suite. Peak list files were searched by the ANDROMEDA a search engine, incorporated into the MaxQuant framework [58], against the IPI-mouse (version 3.68) containing both forward and reversed protein sequences. Initial maximum precursor and fragment mass deviations were set to 7 ppm and 0.5 Da, respectively, but MaxQuant achieved sub-ppm mass accuracy for the majority of peptide precursors. The search included variable modifications for oxidation of methionine, protein N-terminal acetylation and carbamidomethylation as fixed modification. Peptides with at least six amino acids were considered for identification specifying as enzyme LysC allowing N-terminal cleavage to proline. The false discovery rate, determined by searching a reverse database, was set at 0.01 for both peptides and proteins. Identification across different replicates and adjacent fractions was achieved by enabling matching between runs option in MaxQuant within a time window of 2 minutes. Quantification of SILAC pairs was performed by MaxQuant with standard settings using a minimum ratio count of 2.

All bioinformatic analyses were performed with the Perseus software (http://www.perseus-framework.org/↗). To determine the subset of cycling proteins, each protein expression profile is fitted to a cosine with a fixed period of 23.6 h and the amplitude and phase as free parameters. Profiles are ranked by their variance ratio. This is the part of the variance explained by the fit divided by the contribution to the variance that is not accounted for by the fit. Based on this ranking we determine a permutation-based false discovery rate by repeating the same procedure 1,000 times on the same profiles but with scrambled time labels, except for the technical replicates that were preserved. This permutation based procedure is similar to the one applied to FDR calculations for differential expression analysis in [59]. Hierarchical clustering was done in a phase-preserving way by restricting the order of elements to the one determined by the output of the cosine model-based fitting. During the growth of the tree in hierarchical clustering, only those links were permitted that conserve this order. While obviously the order of the terminal branches is not an outcome of the algorithm, the cluster structure is still a non-trivial result of the clustering.

A standard PCA analysis was performed in the Perseus software. The expression data matrix has protein groups as rows and samples as columns and contains logarithms of ratios. Missing values were imputed by drawing random numbers from a normal distribution to simulate signals from low abundant proteins. The width parameter of this normal distribution was chosen as 0.3 of the standard deviation of all measured values and the center was shifted towards low abundance by 1.8 times this standard deviation. These parameters were empirically determined to result in good performance over many different proteomics data sets. The row means were subtracted from the matrix. Then the PCA was performed by singular value decomposition.

Total cell protein extracts were prepared as mentioned above and 50 µg of protein was used for western blotting performed according to standard protocols. Antibodies used were PER21A from Alpha Diagnostics International, RAB1 from Sigma-Aldrich, SARA1; GAPDH and RAB10 from Cell Signaling. The density of the bands obtained in the western blots was calculated with the freely available gel analyzer program ImageJ.

RNA extraction from liver was done using the RNeasy kit according to manufacturer's protocol (Qiagen). cDNA was synthesized from 2 µg of liver total RNA using First Strand cDNA Synthesis kit following the supplier's instruction (Fermentas). Quantitative PCR reactions were done by amplifying ten per cent of the cDNA with Sybr Green master mix (Applied Bioscience) on a CFX96 Real Time System (BioRad). Mean values were calculated from triplicate PCR assays for each sample and normalized to those obtained for Gadph transcript.

Interaction network analysis of the cycling proteome was performed with the STRING search tool (version 9.0) using medium to high confidence (0.5–0.7) and with co-expression and experiments as active prediction methods. Using these parameters we obtained interaction scores for approximately 70% of the liver rhythmic proteins. Data visualization was done with Cytoscape 2.8.2 where we then combined interaction scores with phases of protein abundance.

Data availability: The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium ( http://proteomecentral.proteomexchange.org↗ ) via the PRIDE partner repository [60] with the dataset identifier PXD000601.