L. lactis was grown in continuous culture at different growth rates in the conditions previously described [21] and samples were taken for both transcriptome and proteome analysis at three dilution rates, i.e. 0.09, 0.24 and 0.47 h⁻¹: the lowest growth rate (μ = 0.09 h⁻¹) was chosen as reference. Despite the small size of the L. lactis genome (2310 genes [22]), a total of 346 different proteins were quantified corresponding to 308 different proteins measured in each repetition for each of the 3 steady states. Among these proteins, 193 showed differential profiles in response to a growth rate increase: 88 with reduced level and 105 with higher level. All the proteins displaying a significant level of modification for at least one of the dilution rates are listed in Table 1. In accordance with what has previously been found with transcriptome analysis [21], increased levels of proteins related to biogenesis were observed when the growth rate was increased, i.e. proteins related to transcription (GreA, NusA, QueA, RpoA), translation and more specifically ribosomal proteins (GatA, GatB, RplE, RplI, RplJ, RplK, RplM, RplN, RpmE, RpsA, RpsF, RpsT), enzymes related to fatty acid and phospholipid metabolism (AccA, AccD, FabD, FabF, FabG1, FabH, FabZ1, HmcM, ThiL, YdiD, YscE), two proteins involved in cell division (FtsY, FtsZ), and some proteins associated with purine, pyrimidine, nucleoside and nucleotide metabolism (Add, Adk, Apt, DeoB, GuaA, GuaC, Hpt, NrdE, PydA, PyrC, PyrE, RmlA, RmlB, Upp).

Proteome profiles differed between the various stress-related proteins. On one hand, the two chaperones DnaK and GroEL, the superoxide dismutase associated to oxygen stress SodA, and DpsA, were found in higher quantity, while on the other hand, the cold shock associated protein CspE, ClpC and the adaptation related peroxidase Tpx, had decreased levels in response to growth rate increase. Besides those opposite punctual regulations, other proteins encoding important functions involved in stress protection such as ATPases or peptidases (excepting PepP), were present at constant levels, independently of the growth rate. This lack of general tendency observed here at proteomic level was also observed at transcriptomic level [21]. In contrast, a wide down-regulation of genes involved in stress protection was observed in yeast when growth rate was increased [23],[24]. Finally, one can notice that the two single phage-related proteins measured in those proteomics experiments showed significantly reduced levels at high growth rate. This last observation can be connected with the previously described massive down regulation of the expression of phage-related genes [21].

Transcriptomic and proteomic analyses were performed with cells collected simultaneously from the same fermentor; thus data can be strictly compared. Proteins and their corresponding transcript levels were compared individually. Transcriptomic data were already available [21] but were nevertheless re-processed so as to obtain concentrations rather than abundances (see Materials and methods). For proteomic data, concentration and abundance values are expected to be similar (see Materials and methods). For each growth rate, transcriptomic and proteomic mean values with their standard deviations are given in supplementary data (Table S1). The mRNA/protein ratios were not constant for the different genes since, at a given growth rate, data were spanned among five orders of magnitude (Figure 1C). These variations were linked both to protein and mRNA changes though protein concentrations were globally more spanned than mRNA concentrations (4 and 2 log of magnitude respectively; see figure 1A and 1B). mRNA/protein ratios were compared between two conditions using the lowest growth rate (0.09 h⁻¹) as reference (Figure 1D) and they globally increased with the growth rate. This tendency was confirmed when we analyzed similarly data corresponding to the maximum growth rate of 0.88 h⁻¹. These last data, also available in our group, were obtained in batch culture during the exponential growth phase, since this high growth rate could not be reached in continuous culture without any wash out of the cells from the chemostat [25].

What normally occurs in bacterial cells is the transcription of genes into mRNA, which are then translated into proteins that can be either diluted by growth or degraded (Figure 2). Hence, protein concentration is determined not only by translation rates but also by dilution and degradation, therefore the following balance equation can be written:(1)

Rates of mRNA translation or protein degradation/dilution are assumed to be not constant and related to mRNA or protein concentration respectively. Biological rates were expressed as first order kinetics of their substrate concentration as previously postulated in L. lactis [26],[27] but also in yeast strains [19],[28]. Such modeling approach at the genomic scale is rare in the literature and dynamic experimental data allowing more elaborated kinetics to be hypothesized are not available. Making more complex those rate expressions would thus not make sense today. Dilution constant corresponds to the growth rate (μ), degradation constant (k″) is proportional to protein half-life (t_1/2 = ln(2)/k″) and k′ represents the translation efficiency. At steady state, the various concentrations are expected to remain constant, the time derivative of the protein concentration is equal to zero and the previous equation can be simplified and reorganized as follows:(2)

In the chemostat cultures at the various growth rates, cells are at steady state; similarly, during the exponential growth phase, cells are physiologically stable and are also considered to be at steady-state [29]. The previous observation, establishing a relationship between mRNA/protein ratios and the growth rate for these four steady states (see above), is in accordance with this last equation (2). 171 different mRNA and protein couples were available in each repetition of the various steady states (intersection of 308 couples in the 3 chemostat steady states and 191 in the batch). For only a few proteins were probes missing on the microarray; hence it was not possible to rebuild these couples. In order to estimate translation efficiency and protein degradation rate, the best mathematical solution to the equation (2) was sought, using numerical estimations performed on Matlab. The k′ and k″ values were postulated to be positive, in accordance with biological reality. Different solutions with k′ and/or k″ constant, directly or inversely proportional to μ were investigated. Estimation of the best fitting solution was based on the least square criterion [30]. For the 171 couples, the mean sum of the squared residuals (difference between a ratio and its estimation) associated to every combination are given in Table 2. Considering the lower mean sum of the squared residuals, the best solution was obtained when both k′ and k″ were proportional to 1/μ (k′ = α/μ and k″ = β/μ). Hence the equation (2) could be written as follows:(3)

The mRNA/protein ratios were thus linked to the growth rate (μ) through a polynomial function of order two (μ²), which is consistent with the visual observation of the various curves (not shown). For each mRNA–protein couple the reliability of the two estimated constants α and β was evaluated by their associated R². All regression coefficients are listed in Table 3. The mean linear coefficient (R²) associated to this model was 0.83±0.04. Finally, the consistency of our modeling approach was checked when removing the data of the batch exponential growth phase from the analysis. A high mean R² of 0.77±0.12 was still obtained using chemostat data exclusively. On the contrary, data not at steady state coming from other growth phases in batch cultivation could not be included. Indeed when taking into account mRNA/protein values during growth deceleration or during stationary phase, R² was strongly affected and dropped to 0.21±0.03 and 0.24±0.04 respectively.

The model (1) states that translation rate is proportional to the concentration of mRNA species which assumes that translation is mRNA-limited. An alternative hypothesis, would be the saturation of the ribosome with mRNA, as previously postulated in E coli [31], which implicitly suggests competition of any specific mRNA with all others to be the determining factor in synthesis of the corresponding protein. We have thus tested the model (2) with mRNA abundances rather than concentration values. The modeling approach was robust since similar k′ and k″ dependencies to the growth rate (k′ = α/μ and k″ = β/μ) were obtained (Table 2; data in italic). However the expression of individual protein/mRNA as a function of μ² (3) had generally lower R² (only 48 couples with R²≥0.90 compared to 130 with concentrations), indicating the modeling approach with abundances was less satisfactory than with concentrations.

In order to carry on our modeling approach, the data corresponding to mRNA concentrations were filtered and only couples with R²≥0.90 (130) were retained for further analyses (Table 3). It could be noticed that among the 41 eliminated couples, 15 displayed non monotonous evolutions of their mRNA/protein ratios against μ² and 20 had very low mRNA/protein ratios, which were thus more sensitive to errors.

The k′ and k″ coefficients were numerically calculated for each protein and in each growth condition from the values of α and β (Table 2) by the relation k′ = α/μ and k″ = β/μ. Value distributions (Figure 3) demonstrate the wide variability of k′ and k″ among proteins but also between growth conditions. The k″ decreases when growth rate increases which is consistent with the general idea that protein degradation is high in stationary phases [32]. Protein half-lives were calculated and median values of t_1/2 were respectively 23, 61, 119 and 224 min for 0.09, 0.24, 0.47 and 0.88 h⁻¹. These values are in the same order of magnitude as those obtained recently for S. cerevisiae (mean 43 min at a growth rate of 0.1 h⁻¹ [19]). Like k″, the translation efficiency k′ is also expressed as 1/μ function, indicating that, when the degradation process increases at low growth rate, the translation efficiency is also increasing in order to attenuate this negative biological effect.

Due to the restricted size of the dataset but also to the non-uniform distribution of detected proteins in the various functional categories, it was not possible to use statistical tests to rigorously determine functional enrichments in extreme values of k′ or k″. However, among the 15 genes that are translated the most efficiently (highest α values in Table 3), one can notice the over-representation of genes involved in major cellular processes: the Tig chaperone [33] and proteins involved in replication (HslA, which can unwind DNA and plays a role in its supercoiling, [34]), and translation (ribosomal proteins: RplA, RplF, RplK, RplN, RpsT). Carbon metabolism is also represented by 7 proteins (GapB, EnoA, FbaA, Pmg, Pyk, TpiA, Ldh), all belonging to glycolysis which is the major metabolic pathway for energy production in L. lactis. The extremely stable proteins correspond to null values of β, and consequently k″, were represented by a group of 26 proteins. Remarkably, half of them were related to stress responses: ClpE protease, PepC and PepP peptidases, three reductases that are usually linked to oxidative stress (AhpC, TrxB1, YpjH), but also MurF, involved in parietal structure, YtgH, which is homolog to Staphylococcus aureus alkaline stress protein [35], YtaA and YahB two hypothetical protein sharing homologies with E. coli universal stress protein Usp [36], YuhE, whose E. coli homologue is involved in copper resistance [37], and two cysteine desulfurases (YeiG and YseF) whose corresponding genes in E. coli are involved in oxygen and copper stress responses [38]. Moreover, those extremely stable proteins are rather in the last third for translation efficiency. Thus L. lactis may limit degradation of stress-related proteins so as to maintain a minimal pool ready to use in case of emergency, which is biologically relevant.

Biological determinants of translation efficiency and protein stability were investigated through correlation studies. Correlations providing a Spearman coefficient (R_Spearman) with associated p-values lower than 0.05 were considered as significant. The codon adaptation index (CAI) positively correlates with k′ (R_Spearman = 0.57). Since CAI directly reflects translation efficiency during the elongation step [39], this result validates our translation efficiency estimations. Translation efficiency is also tightly related to the amino acid composition of proteins. A negative correlation of k′ was obtained with tyrosine, cysteine, histidine, aspartic acid and isoleucine frequencies while lysine and alanine richness had a positive influence (Table 4). The amino acids the most used have a positive influence on k′ whereas those with a negative effect are the less frequent ones (Table 4). The single exception is for isoleucine, but since it is the limiting nutrient it is not surprising to find it negatively correlated with translation efficiency, despite its high frequency in L. lactis proteins. This amino acid bias, together with the codon bias (CAI), shows that translation efficiency is strongly dependent of the gene sequence. This optimized functioning state is probably the result of a long evolutionary process. Finally it was found that translation efficiency is affected by protein length: the longer the protein, the more k′ decreases (R_Pearson of −0.18). This negative correlation with length has already been reported for yeast [19] and can possibly be explained by a decrease of the ribosome density on long mRNA as previously shown for S. cerevisiae [40]. The only apparent correlation emerging for protein degradation constants k″ is a negative influence of cysteine richness (Table 4).

Degradation and dilution by growth are both involved in protein disappearance and are competitive reactions. The degradation and dilution constants, k″ and μ, can be directly compared. The k″ is higher than μ at low growth rate but becomes lower after a critical value of 0.39 h⁻¹ (Figure 4). The role of the degradation may thus be major at low growth rate while dilution may become the main phenomenon at fast growth.

More generally, variations in protein concentration between two conditions can be related to changes in the three rates: protein synthesis, degradation and dilution. In order to better understand this regulatory node and identify which are the major controls, the quantitative involvement of the different actors was analyzed. Regulation coefficients corresponding to the protein level control were estimated with a method based on the one developed on S. cerevisiae [41],[42]. Derivation of equation (2) leads to the following relationship:(4)

The term of the equation (4) represents translation control on protein concentration and is called ρ_t while the term −, named ρ_d, includes both the dilution and the degradation and represents protein control by disappearance. ρ_t and ρ_d were estimated for each growth rate interval (between 0.09 and 0.24 h⁻¹; between 0.24 and 0.47 h⁻¹; between 0.47 and 0.88 h⁻¹). The values of ρ_t were used to elucidate the nature of the control and are given in supplementary data (Table S2). If ρ_t≤0, protein disappearance is the major controlling mechanism; if ρ_t≥·1, it is translation; and if 0·< ρ_t·< ·1, the control of protein is shared. The nature of the control for a given protein and its strength differed in the various the growth rate intervals. However a constant control by disappearance was observed for 6 proteins distributed all over the metabolism (Als, GreA, LplL, PyrC RplQ, ThrC, see Table 3 for associated functions). Inversely the unknown protein YpdC was the single one constantly controlled by translation process. Independently of the growth rate, protein levels are mostly controlled by disappearance (Table 5). Translation control strongly decreases, and protein level control becomes less specific and more and more shared with increasing growth rates. Similar conclusions were valid when ρ_d was used for the control analysis instead of ρ_t (data not shown).

Lactococcus lactis ssp. lactis IL1403, whose genome has been entirely sequenced [22], was grown as previously described [21]. Briefly, three different growth rates have been studied, namely 0.09, 0.24 and 0.47 h⁻¹ during anaerobic chemostat cultures (under nitrogen atmosphere and regulated pH) on a chemically defined medium limited by isoleucine concentration. For each steady-state, samples have been harvested in at least quadruplicate with a minimum delay of five doubling time between each sampling.

Transcriptomic data (geo platforms GSE10256 [21] for chemostat culture and GSE12962 for batch exponential phase) were already available. Briefly, these transcriptomic analyses had been obtained with a constant amount (10 µg) of total RNA (mRNA, ribosomal RNA and transfer RNA) labeled by retro-transcription (³³P) and hybridized on nylon membrane as previously described [49]. Three independent biological repetitions were used. These transcriptomic data had been normalized by all spots' mean intensity and thus corresponded to mRNA abundances. They were reprocessed here in order to calculate mRNA concentrations with the method previously described [49]. Raw data were first standardized by the all spots' mean intensities of the reference membrane (and not with its proper membrane) in order to eliminate the bias of the radioactivity level between the various repetitions and then corrected by total RNA concentration in order to take into account changes in intracellular RNA yield in the cells. Consistent with previous results [50], this yield increased significantly with the growth rate in L. lactis (3.58±0.39, 4.92±0.52, 7.34±0.28 and 11.06±0.23 g for 100 g cell dry weight at μ = 0.09, 0.24, 0.47 and 0.88 h⁻¹ respectively). Since the amount of RNA to perform transcriptomic analysis is maintained constant in order to avoid retro-transcription labelling bias, these RNA yield changes are completely hidden by the technology.

The total raw intensity of the membrane without any normalisation represents the amount of mRNA in the RNA sample used for transcriptomic analysis (10 µg). This total intensity was constant at each growth rate and lower than the saturation threshold (mean value of 1660±584, 1462±383, 1389±366, 1474±367, 1496±425 at μ = 0.09, 0.24, 0.47 and 0.88 h⁻¹ respectively). Thus, it can be deduced that the ratio mRNA/total RNA was constant and assuming that ribosomal RNA is the major component of total RNA we can postulate that the fraction mRNA/ribosomal RNA is independent of the growth rate.

For each condition, three repetitions were performed with independent cultures, extractions and electrophoresis. Bacteria were harvested from the cultures and cell pellets were washed twice with ice-cold 200 mM Na-phosphate, pH 6.4 and re-suspended in 4 ml of 20 mM Na-phosphate buffer, pH 6.4, 1 mM EDTA, 10 mM tributylphosphine, a cocktail of protease inhibitors (P8465; Sigma Aldrich, St Louis, MO) 20-fold diluted and catalase 40 U/ml (C3155; Sigma Aldrich, St Louis, MO) to limit isoform formation. The cell suspension (approximately 35 units of optical density at 600 nm [OD600]/mL,) was transferred to the pre-cooled chamber of a BASIC Z cell disrupter (Celld, Warwickshire, United Kingdom) and was subjected to a pressure of 2,500 bars. The suspension was centrifuged at 5,000×g for 20 min at 4°C to remove unbroken cells and large cellular debris. The supernatant was collected and centrifuged at 220,000×g for 30 min at 4°C. The total protein concentration in the resulting supernatant (cytosolic fraction) was determined with the Coomassie protein assay reagent (Pierce, Rockford, IL) using bovine serum albumin as standard and was included between 1 and 2 mg/mL. The cytosolic fraction was aliquoted and stored frozen at −20°C.

R² calculations and equations resolution were perform with MATLAB software.

Correlations were estimated using R free statistical software to calculate Spearman rank correlation coefficient and the associated p-value.