Synechocystis sp. PCC 6803 was cultivated in 1.6 L BG-11 (pH = 7.8) at 30°C in a cylindrical photobioreactor (D = 10 cm, V = 2 L, baffled). The culture was illuminated with an LED light jacket covering the sides of the cylinder (90% red light, 10% blue light). CO₂-enriched air was sparged into the culture (7% CO₂, 330 mL min^–1) and the impeller stirring rate was set to 150 rpm. Cells were grown in turbidostat mode (OD₇₃₀ set point: 0.65–0.80) under an artificial day-night light regime (Day: sinusoidal, max 500 μmol photons m^–2 s^–1; Night: dark) for seven days, at which point the diurnal pattern of dissolved oxygen (growth rate proxy) became stable over subsequent days. Five time points (1 h before/after sunrise, midday, 1 h before/after sunset) were then sampled at -1, 1, 6, 11, 13, 30, 35, 37, 47, and 49 h relative to the first subjected sunrise. Two replicate cultivations were conducted. In the first cultivation, two replicate samples were collected at all five diurnal time points for mRNA sequencing and ribosome profiling. In the second cultivation, one and two replicate samples were collected for ribosome profiling and quantitative proteomics, respectively. Eleven out of fifteen collected ribosome profiling samples were analyzed, which resulted in two replicate measurements at all time points, except 1 h after subjected sunrise which had three. The correlation between ribosome profiling replicates within the same cultivation was similar to the correlation between replicates of different cultivations (r ≈ 0.99 and r ≈ 0.97, respectively), which indicated that results were reproducible across cultivations (Supplementary Figure 1). Cultivation data is shown in Supplementary Figure 2.

Analysis was performed using R v.3.6 scripts¹. Specific growth rates (μ) were determined over the time course of the first cultivation at 30-second intervals according to a mass balance-derived equation of biomass in the culture:

Where D is the dilution rate (h^–1) and OD₇₃₀ is the optical density (cell mass concentration proxy) of the culture. Dilution rates were calculated by dividing time-specific medium feed rates (L h^–1, automatically regulated) with the volume of the culture (1.6 L). To remove noise, each time point was assigned with the 40% truncated average of the closest 480 time points (40% most extreme values removed in each window). The noise filtering step was repeated once. The optical density was measured automatically at 880 nm and converted to OD₇₃₀ using a conversion factor based on offline OD₇₃₀ measurements during the time course. Noise was removed from OD₇₃₀ values in two subsequent filtration steps. In the first step, each time point was assigned with the 40% truncated average of the 960 closest time points. In the second step, each time point was assigned with the average of the closest 360 time points. The change in OD₇₃₀ over time (d⁢O⁢D730d⁢t) was calculated at each time point as the slope in noise-removed OD₇₃₀ (1 h centered time intervals). Noise was finally removed from calculated growth rates by assigning each time point with the 40% truncated average of the 720 closest time points, and negative values were replaced by zero. The average diurnal growth rate was based on growth rates determined across the two days of sampling (0–48 h). Growth rates determined in the second cultivation experiment were prone to error and therefore not reported (higher noise levels and longer data acquisition intervals). Data is shown in Supplementary Figure 2.

Culture medium was removed by centrifugation and cell pellets snap-frozen in liquid nitrogen (stored at –80°C). Cell lysis was performed using lysozyme treatment and vortexing with glass beads. Total RNA was extracted from the cleared lysate with hot phenol/chloroform and isopropanol precipitation and remaining DNA was removed using DNase I. The amount of rRNA was subsequently reduced using the RiboMinus Kit, Bacteria (ThermoFisher, K155004) according to the manufacturer’s instructions. Sequencing libraries were prepared using the NEBNext Ultra Directional RNA Library Prep Kit (New England Biolabs, E7420). Libraries were sequenced on an Illumina NextSeq500 platform (75 bp read length, single end). For details, see Karlsen et al. (2018). Raw sequencing data are available at the European Nucleotide Archive under accession number PRJEB42778.

Cells were rapidly harvested by vacuum filtration and snap-frozen in liquid nitrogen (stored at –80°C). Frozen cells were lysed using cryogenic grinding. The frozen lysate was thawed and cell debris was removed by centrifugation. Polysomes in the lysate were immediately digested with micrococcal nuclease, and generated monosomes were extracted by sucrose gradient ultracentrifugation and fractionation. Total RNA was extracted from monosomes with hot phenol/chloroform and isopropanol precipitation. Ribosome protected mRNA fragments were then extracted by size selection on a denaturing polyacrylamide gel (20–40 nt) and subsequently converted into a sequencing library using the NEBNext small RNA library prep set (New England Biolabs; E7300). Libraries were sequenced on an Illumina NextSeq500 platform (75 bp read length, single end). For details, see Karlsen et al. (2018). Raw sequencing data are available at the European Nucleotide Archive under accession number PRJEB42778.

Culture medium was removed by centrifugation and cell pellets snap-frozen in liquid nitrogen (stored at –80°C). Thawed cell pellets were suspended in 125 μL solubilization buffer (200 mM TEAB, 8 M Urea, protease inhibitor). 100 μL glass beads (100 μm diameter) were added to the cell suspension and cells were lysed by bead beating in a Qiagen TissueLyzer II (5 min, f = 30/s, precooled cassettes). Cell debris was removed by centrifugation at 14,000 × g, 30 min, 4°C, and supernatant was transferred to a new tube. Protein concentration was determined using the Bradford assay (Bio-Rad). For reduction and alkylation of proteins, 2.5 μL 200 mM DTT (5 mM final) and 5 μL 200 mM CAA (10 mM final) were added, respectively, and samples incubated for 60 min at RT in the dark. Samples were diluted 8-fold with 700 μL 200 μM TEAB. For digestion, Lys-C was added in a ratio of 1:75 w/w to protein concentration, and samples were incubated at 37°C and 600 RPM for 12 h. Trypsin was added (1:75 w/w) and samples incubated for 24 h at the same conditions. Samples were acidified with 100 μL 10% formic acid (FA) and insoluble compounds were removed by centrifugation (14,000 × g, 15 min, RT). Peptide samples were then cleaned up using a solid phase extraction protocol (Sep-Pak 1cc 50 mg A C18 cartridges, Waters) according to the manufacturer’s recommendations. Briefly, Sep-Pak columns were equilibrated with 1 mL acetonitrile (ACN) and 2 × 1 mL 0.6% acetic acid. Samples were loaded on columns and washed twice with 1 mL 0.6% acetic acid. Peptides were eluted from the column in 500 μL elution buffer (0.6% acetic acid, 80% ACN) and dried in a speedvac for 2 h, 37°C. Dried peptides were frozen at –80°C and dissolved in 10% FA to a final concentration of 1 μg/μL before MS measurement.

Lysates were analyzed using a Thermo Fisher Q Exactive HF mass spectrometer (MS) coupled to a Dionex UltiMate 3000 UHPLC system (Thermo Fisher). The UHPLC was equipped with a trap column (Acclaim PepMap 100, 75 μm × 2 cm, C18, P/N 164535, Thermo Fisher Scientific) and a 50 cm analytical column (Acclaim PepMap 100, 75 μm × 50 cm, C18, P/N ES803, Thermo Fisher Scientific). The injection volume was 2 μL out of 18 μL in which the samples were dissolved in the autosampler. Chromatography was performed using solvent A (3% ACN, 0.1% FA) and solvent B (95% ACN, 0.1% FA) as the mobile phases. The peptides were eluted from the UHPLC system over a 90 min gradient at a flow rate of 250 nL/min with the following mobile phase gradient: 2% solvent B for 4 min, 2–4% solvent B for 1 min, 4–45% solvent B for 90 min, 45–99% solvent B for 3 min, 99% solvent B for 10 min and 99–2% solvent B for 1 min following re-equilibration of the column at 2% solvent B for 6 min. The MS was operated in a data-dependent acquisition mode with a Top 8 method. The MS was configured to perform a survey scan from 300 to 2,000 m/z with resolution of 120,000, AGC target of 1 × 10⁶, maximum IT of 250 ms and 8 subsequent MS/MS scans at 30,000 resolution with isolation window of 2.0 m/z, AGC target of 2 × 10⁵, maximum IT 150 ms and dynamic exclusion set to 20 s. LC-MS shotgun proteomics data are available at the PRIDE Archive² under accession number PXD023812.

The protein content was quantified in the cell extracts used for LC-MS-MS (Bradford assay). Measured concentrations were normalized to the sample’s cell mass concentration (based on external OD₇₃₀ measurements).

Analysis of sequencing data was conducted using python v.2.7 scripts adapted from Becker et al. (2013), R v.3.4 scripts, and bash commands parallelized using GNU Parallel v.20161222 (Tange, 2011)³. FastQC was used to assess the quality and general features of sequencing datasets (Andrews, 2010). Adapter sequences were trimmed off using Cutadapt v1.18 (Martin, 2011). Base calls with a Sanger quality score lower than 20 were trimmed off the ends of mRNA sequencing reads using Sickle (Joshi and Fass, 2011). Ribosome profiling reads with an average Sanger quality score lower than 25 were removed using Seqmagick v0.6.2⁴. Reads shorter than 6 nt were discarded. Reads that mapped to tRNA and rRNA genome sequences were subsequently removed using Bowtie v.1.2.2 (Langmead et al., 2009). Bowtie was used to map remaining reads to the genome, including plasmids (NC_000911.1 + NC_005229.1 + NC_005230.1 + NC_005231.1 + NC_005232.1 + NC_020289.1 + NC_020290.1 + NC_020298.1). A maximum of two alignment mismatches were allowed. If a read mapped to several locations, only the one best alignment was kept. The read was discarded if it could not be mapped to a unique location in this way. The total number of mapped non-tRNA/rRNA reads was ∼2 million and 34–78 million in mRNA sequencing and ribosome profiling samples, respectively. For each mapped mRNA sequencing read, a read count equal to 1 was distributed evenly over all its aligned genome positions. In contrast, the read count of each mapped ribosome profiling read was assigned to a single genome position, 12 nt upstream of the aligned 3’ end. This assigns the read count to the genome position covered by the A-site of the ribosome (Karlsen et al., 2018). As only ribosome profiling reads longer than 24 nt were counted, the total number of counted mapped reads per sample was between 17 and 68 million. The mRNA/ribosome abundance of a gene (RPKM) was quantified by dividing the read count on the gene’s coding sequence with the length of the coding sequence (in 1,000 base pairs) and the total number of counted reads on all coding sequences (in million). Coding sequences were defined according to GenBank files for the NCBI reference sequences NC_000911.1, NC_005229.1, NC_005230.1, NC_005231.1, NC_005232.1, NC_020289.1, NC_020290.1, and NC_020298.1.

Thermo raw spectra files were converted to the mzML standard using Proteowizard’s MSConvert tool (Adusumilli and Mallick, 2017). Peptide identification and label-free quantification were performed using OpenMS 2.4.0 in KNIME (Röst et al., 2016). The KNIME pipeline for MS data processing was deposited on https://github.com/m-jahn/openMS-workflows (labelfree_MSGFplus_Percolator_FFI.knwf). MS/MS spectra were subjected to sequence database searching using the OpenMS implementation of MS-GF+ (Granholm et al., 2014) with the Synechocystis sp. PCC 6803 reference proteome as database (as of 04 April 2019). Carbamidomethylation was considered as a fixed modification on cysteine and oxidation as a variable modification on methionine. The precursor ion mass window tolerance was set to 10 ppm. The PeptideIndexer module was used to annotate peptide hits with their corresponding target or decoy status, PSMFeatureExtractor was used to annotate additional characteristics to features, PercolatorAdapter was used to estimate the false discovery rate (FDR), and IDFilter was used to keep only peptides with q-values lower than 0.01 (1% FDR). The quantification pipeline is based on the FeatureFinderIdentification workflow allowing feature propagation between different runs (Weisser and Choudhary, 2017). MzML files were retention time corrected using MapRTTransformer, and identifications (idXML files) were combined using the IDMerger module. FeatureFinderIdentification was then used to generate featureXML files based on all identifications combined from different runs. Individual feature maps were combined to a consensus feature map using FeatureLinkerUnlabelledKD, and global intensity was normalized using ConsensusMapNormalizer (by median). Protein quantity was determined by summing up the intensities of all unique peptides per protein using ProteinQuantifier.

Analysis was performed using R v.3.6 scripts¹. Genes with less than 30 and 60 reads were initially removed from the mRNA sequencing dataset and ribosome profiling dataset, respectively. In the proteomics dataset, inaccurately measured proteins were removed by discarding those with a log₂ intensity standard deviation greater than 1 at any time point. Only genes with at least two replicate measurements across all time points and across all three datasets were analyzed. The total number of genes (n = 1,126) was mainly limited by the proteomics dataset. Abundance values for each gene were log₂ transformed and then centered around the gene’s daily average log₂ abundance to reflect relative fold changes. To identify genes with diurnal changes in mRNA abundance (considered “cyclic”), differential abundance between time points was analyzed with one-way ANOVA. A gene’s mRNA abundance was considered to change significantly over the day-night cycle if (1) the Benjamini-Hochberg adjusted q-value was less than 0.1, and if (2) the absolute log₂ fold change was greater than 1, between any two time points. Cyclic genes were then clustered into four groups (G1–G4) according to their diurnal mRNA abundance pattern using hierarchical clustering (R function: hcluster; distance measure: Pearson correlation; linkage method: Ward). Maximum cluster separation was obtained when choosing a cluster number of 2 and 4 (average silhouette width of 0.63 and 0.49, respectively). Four clusters were chosen as relatively unique diurnal patterns were visible in each, despite lower average silhouette width. Non-cyclic genes were assigned to a fifth group (G0). The peak-to-peak relative amplitude in mRNA, ribosome and protein abundance was calculated as the maximum log₂ fold change between sample mean values across the day-night cycle. The median relative amplitude of cyclic genes in each dataset was used to summarize and compare the overall relative amplitude observed at transcriptional, translational and protein level. Differential log₂ protein abundance across time points was assessed for cyclic and non-cyclic genes using one-way ANOVA. A Benjamini-Hochberg adjusted q-value less than 0.1 was considered significant.

The change over time of an arbitrary gene’s (J) protein concentration (P_J) was expressed according to the cellular mass balance of that protein:

Where F_{S, J} is the fraction of total bulk protein synthesis (S_TOT) dedicated to protein J, F_{P, J} is the fraction of the total cellular protein concentration (P_TOT) made up by protein J, μ is the growth rate, and k_{D, J} is the gene-specific degradation rate of protein J.

In a similar manner, the rate change of P_TOT was expressed according to the cellular mass balance of total protein:

Where k_{D, MEAN} is the bulk protein degradation rate.

Under the assumption of a constant P_TOT, the term S_TOT is constrained to be proportional to the sum of μ and k_{D, MEAN}:

Substitution of Eq. 3 into Eq. 1 yields:

For the model to better reflect the abundance fraction of protein J measured by protein mass spectrometry, the expression P_J = F_P,J⋅P_TOT as substituted into Eq. 4:

Considering an “average” protein J for which k_{D, J} = k_{D, MEAN}, Eq. 5 can be further simplified to:

Simulations of diurnal protein abundance oscillations were performed using R v.3.6 scripts¹. The model equation (Eq. 5) was solved numerically using an ordinary differential equation solver (R function: ode). F_S,J was set as a function of time with an amplitude fold change corresponding to the experimentally observed median value:

The behavior of F_{P, J} over time in response to different protein turnover scenarios was analyzed by altering the settings of the remaining input parameters μ, k_{D, MEAN} and k_{D, J}. To simulate anticorrelated degradation vs. synthesis, k_{D, J} was expressed as a sine function with a 12 h phase shift relative to F_S,J:

A time depended growth rate was modeled by expressing μ as a sine function during day time:

and as zero during night time.

To investigate whether translational regulation causes protein levels to remain constant during day-night cycles in cyanobacteria, we performed genome-wide measurements of the transcriptome, translatome, and proteome in Synechocystis using mRNA sequencing, ribosome profiling, and quantitative shotgun proteomics, respectively. Cells were adapted to an artificial day-night regime for seven days in a controlled turbidostat culture and samples were collected over separate days at five time points: 1 h before and after artificial sunrise, midday, and 1 h before and after sunset. The maximum and average growth rate was 0.05 and 0.018 h^–1, respectively, and correlated with the light intensity curve (Supplementary Figure 2).

A total of 1126 genes were analyzed, which had at least two replicate measurements across all time points in all three datasets. Of these, 43% showed cyclic diurnal mRNA expression (|log₂FC| > 1, FDR < 0.1, Supplementary Table 1) which is in agreement with microarray-based transcriptomics studies from Synechocystis and Synechococcus elongatus PCC 7942 grown in diurnal light conditions (Guerreiro et al., 2014; Saha et al., 2016). These genes were designated “cyclic genes” and clustered according to their diurnal mRNA abundance pattern into four groups: “G1-G4.” Non-cyclic genes were assigned to a fifth group: “G0” (Figure 1A). Protein synthesis rates, inferred from the number of translating ribosomes, correlated well with cyclic transcription patterns (r = 0.88, Figure 1A). This implies a low degree of translational regulation and confirms that protein synthesis rates oscillate significantly over the day night cycle, in concert with transcript levels. In contrast, oscillation patterns could not be distinguished at the protein level which remained relatively constant. Since the variance between time points was low relative to the variance between replicates, no significant change in protein abundance was found across any time point (ANOVA, FDR < 0.1, Supplementary Table 1). This further explains the low correlation (r = 0.06) between protein synthesis and protein abundance patterns in this dataset (Figure 1A).

To analyze the peak-to-peak oscillation amplitude at the level of transcription, translation and protein, the log₂ fold change between the minimum and maximum time point mean abundance of mRNA, translating ribosome, and protein was calculated for each cyclic gene (termed “relative amplitude”). The median relative amplitude across genes was compared to quantify the overall amplitude reduction from protein synthesis to protein abundance (Figure 1B). The median amplitude was two times lower at level of protein (1.5-fold) compared to the level of transcription and translation (3.0-fold), which is similar to the 2.3-fold median ratio between transcript and protein oscillations reported by Waldbauer et al. (2012). However, the median protein amplitude was probably overestimated here since the error of time point means was high relative to the variation between time point means. For example, if the error of time point means are large and the true protein amplitude is small for a gene, the measured variation across time points will be mostly noise. Consequently, the calculated relative amplitude is likely to be mostly noise, as it will be determined from the maximum fold change across five error-prone time point means. The variance between replicates was larger than the variance between time points for 41% of cyclic genes (ANOVA, SSB/SSW > 1). Thus, the determined median relative protein amplitude is likely to provide a certain over-estimated relative amplitude. Even though gene-specific protein amplitudes could not be determined with precision, there was a trend that proteins with high turnover rate (daily mean protein synthesis rate/daily mean protein abundance) had stronger oscillations (Figure 1C).

In conclusion, our multi-omics analysis shows that the decrease in oscillation amplitude between the mRNA level and the protein level is not caused by translational regulation of protein synthesis. The ratio between the median mRNA oscillation amplitude and the median protein oscillation amplitude was estimated to be more than twofold.

Protein concentrations in cyanobacteria remain largely constant over diurnal cycles, despite significant fluctuations in transcription and protein synthesis. A possible and intuitive explanation is that diurnal peaks in synthesis are counteracted by increased protein degradation. At the same time, this seems unlikely in an evolutionary context, as it implies an ineffective use of cellular resources which would result in decreased fitness. Therefore, we sought to determine whether this observation could solely be the result of a slow protein turnover. For this purpose, we applied a mass balance-based model that describes the change in an arbitrary protein’s concentration in response to diurnal synthesis oscillations (see section “Materials and Methods”). The model takes into account the synthesis, degradation, and growth dilution of the modeled protein (protein J) as well as the synthesis, degradation, and growth dilution of the bulk proteome. The model assumes that the total cellular protein concentration is constant, which has been shown experimentally in cyanobacteria over a range of growth rates and genetic perturbations (Touloupakis et al., 2015; Zheng and O’Shea, 2017). This assumption constrains bulk protein synthesis to be proportional to protein depletion, which is the sum of two processes, bulk protein degradation (described by k_{D, MEAN}) and dilution by cell growth (described by μ). Bulk protein turnover, defined as bulk protein synthesis rate divided by bulk protein abundance, is thus also proportional to protein depletion (see section “Materials and Methods,” Eq. 3). We measured the total protein concentration in cell extracts and found it to be constant across time points (p ≥ 0.4, Supplementary Figure 3). In any case, deviations from this assumption do not have a significant impact on the amplitude of protein oscillations as long as the change is restored within the time span of the day-night cycle (see section “Discussion”).

The effect of protein turnover rate on diurnal protein oscillations was explored using the mathematical model described above. Model parameters were selected to simulate biologically relevant cellular scenarios. The synthesis rate of protein J was set to oscillate with an amplitude equal to the median amplitude determined by ribosome profiling (Figure 1B) and μ was set to the observed daily average (Supplementary Figure 2). The bulk protein degradation rate was set to the median degradation rate reported for microalgae and plants (0.01 h^–1), as it has not been determined experimentally in cyanobacteria (Table 1). In a first simulation, the gene-specific degradation rate of protein J (k_{D, J}) was set equal to the bulk degradation rate so as to mimic the response of an “average gene.” This resulted in a relative protein amplitude of 1.11-fold, which is similar to the relative protein amplitude determined from the experimental data (Figure 2A). Increasing the bulk turnover rate (μ + k_{D, MEAN}) from 0.028 to 0.118 h^–1 (by increasing k_{D, MEAN} from 0.01 to 0.1 h^–1) resulted in increased relative amplitude (1.11–1.53-fold) and decreased lag time of a protein J’s oscillations. The change in relative amplitude was in this case caused by changes in the absolute protein abundance difference between peak and trough (termed “absolute amplitude,” see legend Figure 2). A positive correlation between protein amplitude and bulk protein turnover was also predicted implicitly in model Eq. 6, where a higher bulk turnover increases the protein response dP/dt, which leads to a faster change and increased amplitude (see section “Materials and Methods”). Equation 6 further shows that the direction of the protein change is determined by the difference in protein synthesis fraction and protein abundance fraction (F_{S, J}–F_{P, J}). The abundance fraction will therefore become equal to the synthesis fraction over time, if the synthesis rate of a gene J is constant (e.g., at steady state growth) and k_{D, J} = k_{D, MEAN}. More importantly, this implies that a change in synthesis rate from one steady state to a new one, will result in an abundance change that is at most equal to the synthesis change, if given enough time to reach the new steady state (∼5 protein half-lives). Thus, under a diurnally changing synthesis rate, the protein amplitude is bound to be less than (or at most equal to) the protein synthesis amplitude, unless the protein half-life is much shorter than the time period of the day-night cycle (i.e., relatively high protein turnover).

A positive effect on protein J’s relative amplitude was also observed when its gene-specific protein turnover rate was high relative to the bulk protein turnover rate (by setting k_{D, J} > k_{D, MEAN}, Figure 2B). In this case, the increased relative amplitude was a result of reduced daily mean abundance of the modeled protein (and unchanged absolute amplitude), instead of increased absolute amplitude as in Figure 2A. Similarly, a gene-specific turnover rate lower than the bulk degradation rate (k_{D, J} < k_{D, MEAN}) results in increased daily mean abundance and reduced relative amplitude (data not shown). A positive correlation between the daily mean protein synthesis to abundance ratio and the gene-specific protein turnover, as determined by k_{D, J} in the model, is consistent with the general definition of protein turnover (turnover rate = synthesis rate/protein abundance, at constant protein abundance). Thus, the model predicts a positive trend between the relative protein amplitude and the gene-specific protein turnover rate, which possibly explains the experimentally observed positive trend between relative protein amplitude and gene-specific turnover (Figure 1C). The model further suggests that high-amplitude proteins are likely to possess a high gene-specific degradation rate relative to the bulk protein degradation rate. Increased absolute amplitude was also observed when the gene-specific degradation rate was actively regulated (time-dependent) and anticorrelated to the synthesis rate (Figure 2C). Thus, particularly high relative amplitudes are possible for a subset of genes even at slow bulk protein turnover, if regulated degradation is fast enough to also reduce the daily mean abundance. Allowing the growth rate to fluctuate according to our experimental data altered the pattern of protein oscillations but did not have a significant effect on the protein amplitude (Figure 2D).

These results demonstrate that the relative amplitude of a protein depends on the bulk protein turnover and the protein’s specific turnover. The bulk protein turnover acts on the protein’s absolute amplitude (and all other proteins), while the latter acts on the daily mean abundance of the protein. The model further shows that the observed reduction in oscillation amplitude of a given protein can be attributed solely to a low bulk protein turnover, corresponding to the experimentally determined growth rate of 0.018 h^–1 and a bulk protein degradation rate of 0.01–0.05 h^–1.