Nonlethal deleterious mutation–induced stress accelerates bacterial aging

Significance Genomic heterogeneity is pervasive in living organisms, from bacterial ecosystems to human tissues, allowing both adaptability and/or reestablishment of homeostasis in the proper context. In nature, mutations that cause loss of gene function are a major contributor to bacterial ecosystems’ adaptation to new environments. However, how bacteria handle random, nonlethal, but deleterious loss-of-gene-function mutations is only partially understood. In this study, we show both at population- and single-cell levels that single-gene-deleted Escherichia coli mutants activate a homeostatic stress response, i.e., allostasis, which reflects both the nature of the environment and the deleted gene function. Stress and allostasis, however, also lead to early onset of aging in mutant bacteria, thus representing a functional cost for maintaining population-level resilience to environmental challenges.


Figures S1 to S8
Tables S1 to S2 Legends for Movies S1 to S2 SI References

Other supporting materials for this manuscript include the following:
Movies S1 to S2

Supporting Information Text 1. Environment-dependent growth differences in deletion mutant strains
When grown in LB medium, some mutants with a deleted ATP synthase subunit (ΔatpA, ΔatpD, ΔatpE, ΔatpF, ΔatpH), a DNA repair and stress response related gene (ΔrecB), ubiquinone biosynthesis related genes (ΔubiF, ΔubiH), a gene related to electron transport (ΔnuoJ), a tricarboxylic acid cycle related gene (ΔsucB), a bacteriocin transport gene (ΔtolR), and one other strain (ΔfkpB), were found to have significantly lower growth rates compared to wild-type (wt) E. coli (Figure S2B).However, of these strains, only six of them showed a significant lower saturation point compared to wt cells (ΔatpA, ΔatpD, ΔatpH, ΔrseA, ΔsucB, ΔubiF,).When grown in M9CG or M9G media, a different set of mutants showed significantly lower growth rates.These strains consist of ATP synthase subunit deleted mutants (ΔatpA, ΔatpD, ΔatpE, and ΔatpF), mutants with deleted genes related to sulfur metabolism (ΔcysC, ΔcysD, ΔcysE, ΔcysG, ΔcysH, ΔcysI, ΔcysN, ΔcysP, ΔcysU, ΔcysW), and other mutants (ΔrimM, ΔfkpB, ΔiscS, ΔlipB, Δlpd, ΔnuoA, ΔnuoK, ΔrecB, ΔydjI, ΔubiF, ΔubiH,).Seven of these strains (ΔcysC, ΔcysD, ΔcysG, ΔcysH, ΔcysI, ΔcysN, and ΔcysU), demonstrated a lower saturation point than wt E. coli cells in M9CG.Also, all of these mutants in addition to one more strain (ΔcysE) exhibited a significantly lower SPG when grown in M9G.One strain had lower growth rate than wt in M9CG, but higher in M9G (ΔydjI).As expected, all values for growth rate are lower than wt (except for the two cases of ΔdegP and ΔydjI in M9G), illustrating that no single gene deletion increases the normal growth rate of E. coli cells in these experimental conditions.It is also notable that some mutants have a closer growth rate value to wt E. coli in M9CG than in LB; for example, ΔatpA exhibits a more similar growth to wt E. coli when in M9CG, and even more similar when grown in the nutrient limited medium, M9G.

Gene expression profiles in deletion mutant strains
The E. coli strain BW25113 harbors a complex genetic network comprising of 4,487 genes and 459 pathways 1 , which are responsible for carrying out various biological processes.In this section, we aimed to investigate the response of single gene deleted mutants to stress induced by nutrient limitation.To this end, we analyzed and compared the gene expression profiles of relevant genes of ATP subunit deleted mutants and the wt strain under three different media conditions.Principal component analysis was performed separately within each growth medium, and the first two principal components are presented in Figure S4A.Distinct clustering patterns are evident across the three media conditions, emphasizing significant transcriptional differences between these two strains.
Previous studies have identified genes that become essential when cells were grown in nutrient limited conditions 2 .Figure S4B illustrates the expression profiles for 124 of these genes that were grouped according to their involvement in metabolic pathways 3 .Notably, significant alterations were observed primarily in genes associated with central metabolism, as well as those related to vitamin and cofactor metabolism pathways.These alterations in gene expression likely reflect a cellular response aimed at optimizing resource allocation, rerouting metabolic pathways, and enhancing nutrient uptake and utilization, thus triggering the allostasis response.The identification of these genes provides insight into the regulatory networks that enable cells to maintain essential metabolic functions and survival under challenging nutritional conditions.Furthermore, to investigate how stress response in single cells can lead to different cellular morphologies, we examined the genes related to filamentation and the SOS response in E. coli cells.sulA, minCD, and ftsZ are key genes that play essential roles in the filamentation process 4,5 .Figure S4C displays the expression profiles of these genes, along with other filamentationassociated genes identified as essential 6,7 .The observed downregulation in the expression of ftsZ, which is vital for Z-ring formation and cell division regulation, as well as the upregulation of minE expression, known to restrict Z-ring formation, under nutrient-limited conditions, agrees with the increase in filamentation observed in our single cell experiments.This concurrence between gene expression changes and the alteration in death phenotypes supports the notion that stress conditions can trigger shifts in cell death phenotypes rising from the adaptive mechanisms underlying allostasis.

Probabilities of single cell death and their impact on population growth
Many mathematical models have been proposed for describing bacterial population growth; however, most do not take single cell behavior into consideration 8 .Among these models, empirical analytical equations such as the logistic model are commonly used.In our study, we observed that our population growth curves (from Figure S1 , where  represents the population number as a function of time, ,  0 is the initial cell count,  is the growth rate coefficient, and  is the suppression rate coefficient 9 .Figure S7A demonstrates a good fit of our data to this model.The parameter  is typically associated with growth suppression resulting from nutrient limitation.Here, we aim to clarify whether the growth suppression could be attributed to a combination of growth arrest and death due to bacterial aging rather than due to nutrient limitation. We propose a model that integrates the results of our single cell experiments to predict bacterial death and growth arrest associated with aging.The model is based on two of our key findings: (1) the occurrence of a post-replicative lifespan preceding cell lysis, and (2) the survival functions shown in Figure 5E.The underlying principles of our model are presented in Fig. S7B.Each time-step, denoted as , represents the generation time, which is considered to be constant for simplicity.The fate of each cell is independent of other cells.Within each time-step, a cell faces two possibilities: it can either be in its replicative lifespan and divide with probability   (), or it can enter its post-replicative lifespan with probability   ().In the event of cell division, the same two options are presented in the subsequent time-step (  ( + 1) and   ( + 1)).Conversely, if the cell enters its post-replicative lifespan, it encounters two new options in the next time-step: either it remains in its current state with a probability   ( + 1), or it undergoes cell death with a probability   ( + 1).The values of these probabilities vary across strains and media conditions and are determined based on the survival functions presented in Figure 5E.The red lines in these graphs represent fitted Gompertz distributions, which serve as the basis for acquiring the probability distributions presented in Figure S7C.Let  1 () be the survival function for lysed cells (Figure 5E), and  2 () be the survival function for replicating cells (Figure S6D).By using these functions, we can find the three probabilities as follows:   () =  2 (),   () = 1 −  1 (), and   () =  1 () −  2 ().
Figure S7C presents the outcomes of our model, revealing that bacterial aging has an influence on the population growth experiments conducted in this study.A previous study has investigated single-cell behavior during population growth and introduced three probabilities for individual cells: division, death, or remaining alive but non-dividing 10 .These probabilities are similar to the concepts we have introduced in our model, with the difference that they were attributed to nutrient depletion during population growth.However, other studies have demonstrated that under nutrient-depleted conditions, cells enter a state of growth arrest rather than undergo cell death 11 .
Here, we hypothesize that during population growth, cells undergo these three fates, which occur due to aging.This study reveals that nutrient limitation accelerates the aging process, as evidenced by the survival functions in Figure 5E.Particularly in microbatch experiments, where nutrient depletion occurs rapidly, the aging rate of bacteria can undergo rapid changes due to nutrient scarcity.Hence, the parameter  in the logistic model could be attributed to aging due to nutrient limitation.The precise mechanism by which nutrient depletion affects the aging rates needs further investigation.

Data analysis Microbatch culture data analysis
In a previous study, growth of non-lethal E. coli isogenic single-gene knockout mutants (the Keio collection) on semi solid agar was investigated, and three E. coli colony growth valueslag time of growth (LTG), maximum growth rate (MGR) and saturation point of growth (SPG)were measured using a high-throughput measurement system.The results were able to identify a subset of genes whose deletion delayed colony formation or changed saturation point and presented a better understanding of the range of phenotypes that exist in this strain 12 .We chose 65 mutants among this subset which are listed in Table S1, and measured their MGR and SPG, as explained below.We did not measure LTG for our microbatch culture experiments since the lag time depends on many initial condition factors that are not related to the interests of this study.However, we measured the Area Under Curve (AUC) and time of exponential phase (  ) instead.
To determine the maximum growth rate (MGR) and saturation point of growth (SPG) for each mutant, we employed the following procedure.First, any  600 values below zero were adjusted to zero, as negative values lack meaningful interpretation and result from measurement errors in low concentrations.Subsequently, background measurements obtained from blank samples, where only the growth medium was present in the well, were subtracted from each population growth curve on a point-by-point basis.To identify the MGR, we located the maximum point of the first derivative of  (   0 ), and then selected five points preceding and following this maximum.
A linear regression was fitted to these eleven points, and the slope of the line provided the MGR.For SPG determination, the last 20-50 points of the growth curve were forced to fit a horizontal line, and the value of this line was recorded as the SPG.
To measure AUC and   , we applied methods described in a previous study 13 .We determined the AUC by using MATLAB's trapezoid integration function with time intervals set at 0.1 hours.  was calculated by first applying a piecewise linear fit to  (   0 ) to find the exponential phase.Subsequently, we determined the duration during which the graph exhibited exponential growth.
A demonstration of these calculations for a sample growth curve (wild-type in LB) is presented in Figure S2A.The extracted values from these calculations are presented in Figure S2B-C.

Data reduction and clustering
We calculated the Pearson correlation coefficients between the parameters (Figure S2D), which revealed strong correlations among these variables.This suggests the possibility of reducing the dimensionality of the dataset.To extract the underlying dimensions characterizing the growth curves of microbatch cultures, we performed principal component analysis (PCA).All the parameters (MGRs, SPGs, AUCs and   s) distinguishing the different growth curves of all 65 mutants + 1 wild-type strain across three different media conditions form a 66 × 12 matrix, which we label .Using the built-in MATLAB function "pca", we computed the principal components of this matrix.The scree plot displayed in Figure S2E indicates that more than 95% of the variance can be explained by the first two principal components, justifying a dimension reduction to two.The projection of the entire dataset onto these two principal components is presented in Figure 1A. Figure S2F demonstrates the most relevant parameters to the first two principal components.For clustering analysis, we utilized the hierarchical method and the built-in MATLAB function "clusterdata", with the maximum number of clusters set to two.

Microarray Analysis
Raw Affymetrix CEL files (E. coli Genome 2.0) were processed with Expression Console (version 1.4.1.46),using the Robust Multiarray Average (RMA) to produce normalized probesetlevel expression values and the MicroArray Suite 5.0 (MAS5) algorithm to produce Absent/Present calls for each probeset in each sample.Array quality was also assessed by computing Relative Log Expression (RLE) and Normalized Unscaled Standard Error (NUSE) values using the affy (version 1.70.0) and affyPLM (version 1.68.0)R packages.Principal Component Analysis (PCA) was performed using the prcomp R function with expression values that had been normalized across all samples to a mean of zero and a standard deviation of one.All analyses were performed using the R environment for statistical computing (version 4.1.2).

Single-cell data analysis
The movies acquired from the mother machine experiments were manually analyzed to count the number of dead cells, time of death, percentage of death phenotypes, and number of generations presented in Table S2 and the data points in Figure 5E and S6D.
We note that during the experiment, a subset of cells loses its plasmids and dies as a result of losing the antibiotic resistance gene contained in that plasmid (last column of Table S2).This mode of death was identified by a gradual decrease in the GFP-Fis fluorescence in the cell, while the cell continued to divide (Movie S2).GFP-Fis is produced from the antibiotic resistancecarrying plasmid in the cell.Its sudden loss would lead to a halt in GFP-Fis production and consequently a gradual decrease in the fluorescence intensity in the cell as the cell depletes the protein during subsequent division events.Since the death of such cells is due to our experimental conditions, we did not consider them in our analysis of death phenotypes.
The fits in Figures 5E, S6D, and S8C were performed in MATLAB to the equation  =  (−   (  − 1)) where  is a constant,  is the aging rate, and  is time, using a standard linear least squares regression method.The values for aging rate reported throughout the entire study are the values that were acquired for  using this fit.
The software Oufti 14 , along with custom designed MATLAB codes were used to measure the length of the mother cell.An example output of these measurements along with an explanation of parameter extraction is provided in Figure S6B.The values for exponential elongation rate (EER) and birth length (BL) for each generation were extracted from fitting the data points of each generation to the equation   () =   (0)    , where  is length,  is generation number,  is time, and  is EER.Doubling rate (DR) was determined by reversing the generation time,  = 1   , where   is the generation time.In cases where the individual cell divided one last time and the daughter did not grow in size (which means the last cycle had an exponential elongation rate of zero), generation zero was considered to be the cycle preceding the last division.

Fig. S1. Population growth in microbatch cultures
Growth curves for wild-type E. coli and its 65 tested isogenic mutants in LB (black), M9CG (blue) and M9G (red) growth media are shown.The y-axis is ln ( ), where  is OD600 measurements and  0 is the initial OD600 values.The decline of some growth curves after reaching saturation indicates a death phase of the population.

Fig. S2. Growth parameters of wt and mutant E. coli strains in microbatch cultures
(A) The Maximum Growth Rate (MGR), Saturation Point of Growth (SPG), Area Under Curve (AUC), and time of exponential phase (  ) were extracted from the growth curves using established methods.MGR was calculated as the maximum slope of a fitted line from 5 points before and 5 points after the maximum of the first derivative of (   0 ), where  is  600 measurement and  0 is the initial  600 value.For wt E. coli in LB, the MGR was 0.89 ℎ −1 .SPG was determined by fitting the last 20-50 points of the growth curve to a straight horizontal line.For wt E.coli in LB medium the SPG was 2.65.AUC was determined by integrating the growth curve.For wt E. coli in LB medium the AUC was 32.51.  was found by using a piecewise linear fitting method and finding the time in which the growth is completely exponential.For wt E. coli in LB medium the   was 1.83 ℎ.(B) The MGR, SPG, AUC, and   values that were extracted using the above method for all 65 single-gene mutations presented in Figure S1 are shown here in a color-coded format, colors present the values of the parameters.(C) Plots of MGR against SPG, and AUC against MGR and SPG values reveal high correlations between the variables justifying a reduction in the number of parameters and principal component analysis.(D) Pearson correlation coefficients were calculated for all parameters extracted from the growth curves in three different environments for 65 E. coli single gene deletion mutants and the wt.Some of the correlations are significant justifying a principal component analysis.(E) Cumulative sums of eigenvalues were calculated to determine the percentage of variance explained by each principal component.The first two components explain approximately 95% of the variance of the data.(F) The orthonormal principal component coefficients for each variable and the principal component scores are visualized showing that the most relevant parameters for PC1 are AUC-M9CG and AUC-M9G and for PC2 are AUC-LB,   -M9G and SPG-LB.

Fig. S3. Cysteine metabolism pathway and its enzymes in E. coli
The growth of wt (darker lines) and single-gene subunit deletion mutants (lighter lines) are displayed in graphs surrounding the cysteine metabolism pathway and its enzymes.Seven single deletion mutant strains (∆, ∆, ∆, ∆, ∆, ∆, ∆) belonging to cluster I from Figure 1A, display slow or no growth in nutrient-limited conditions (M9CG and M9G).The other two mutants (∆, ∆) belong to cluster II and display adaptive behaviors in M9CG.The y-axis is ln ( ), where  is OD600 measurements and  0 is the initial OD600 values, and x-axis is time (hours).

Fig. S4. Gene expression profiles across environmental conditions
(A) Principal component analysis was employed to assess the transcriptome dataset under specific media conditions.The values of the first two principal components are presented for each environmental condition, providing an overview of gene expressions in different media.(B) Expression profiles of those 124 genes that become essential when cells are exposed to nutrientlimited conditions and experience stress 2 .These genes are categorized based on their involvement in distinct metabolic pathways, highlighting their functional relevance in response to nutrient limitation-induced stress.(C) Expression profiles of 15 genes essential for filamentation 4- 7 .These expressions align with the observation of an increase in the occurrence of filamentation in the last cell cycle (phenotype II) under nutrient-induced stress.

Fig. S5. Comparison of the growth curves for plasmids-containing and plasmids-lacking strains
(A) A comparison of the growth curves of the two strain types investigated in this study.In the microbatch experiments, cells did not contain any plasmids (w/o plasmid), while in the mother machine experiments, cells contained two plasmids (w/ plasmid), namely pZA3R-mcherry and pRJ2001-GFP-Fis (see Materials and Methods for details).Y-axis is ( ), where  is  600 measurement and  0 is the initial  600 value.(B) A comparison of the MGR, SPG, AUC, and Texp of the two strains.All growth parameters are linearly proportional to each other.The line represents a linear regression of the data points, with the slope value depicted in the plot. .(D) Survival functions of wt and ∆ in three media conditions, where only the aRL of cells are considered as their lifetime.Dotted lines represent experimental data, and red lines are a fit to a Gompretz equation, which allows for estimation of the aging rate ( (ℎ −1 )) shown on each graph.The numbers in the parentheses present 95% confidence bounds., the cell has two options.With probability   ( = 1) it will enter its post-replicative lifespan, and with probability   ( = 1) the cell will divide into two daughter cells.If the mother cell divides, the two daughter cells that are now in generation two (g = 2), will have the same possible fates, but with new probabilities (  ( = 2) and   ( = 2)).If the cell enters its PRL, it will have the probabilities to either die in the next generation (  ( = 2) ) or stay in its post-replicative lifespan (  ( = 2)).(C) Probability functions in all conditions.At the beginning of a cell's lifespan, the cell has a high probability to divide, and a low probability to die or enter its post-replicative lifespan.The probability to enter PRL increases during the middle of its lifespan, and the probability of death increases during the end of its lifespan.5A since some cells were not followed until the end of PRL and could not be included in that dataset.Some experiments did not include the GFP-Fis plasmid, hence it was not possible to differentiate between type Ia and type Ib.The percentage of non-filamented cells are reported as one single number for type I for those experiments.The dashed lines in the last column are experiments where the cells did not contain the GFP-Fis plasmid.For an example of a cell losing its plasmid see Movie S2.Movie S2 (separate file).An example of a cell that loses its PRJ2001-GFP-Fis plasmid and undergoes lysis due to the antibiotics present in the media.This cell loses its plasmid at 54:15 hr, proceeds to divide seven additional times, becomes filamented, and undergoes lysis at 56:15 hr.
) were well-fitted by a logistic model represented by the equation ()

Fig. S6 .
Fig. S6.Microfluidic setup and lifespan distributions of E. coli strains (A) The mother machine was used to track individual bacterial cells for multiple generations until they died.This microfluidic device contains microchannels of similar width of a single wt E. coli cell grown in LB (~1 ) and a length sufficient for multiple cells (~30 ), with a larger channel (width and height of ~30 ) for continuous media flow to support the growth of trapped cells and wash away daughter cells.(B) Example of cell length measurements over generations, with blue dots representing measurements and red lines showing exponential fits to the equation   () =   (0)    , where  is cell length,  is time,  is generation number, and  is the exponential

Fig. S7 .
Fig. S7.Aging probabilities of E. coli determined from single-cell experiments (A) The population growth curves exhibit a good fit to the logistic model.Solid lines represent experimental data and dotted lines depict the logistic model fit.(B) Starting at generation one (g = 1), the cell has two options.With probability   ( = 1) it will enter its post-replicative lifespan, and with probability   ( = 1) the cell will divide into two daughter cells.If the mother cell divides, the two daughter cells that are now in generation two (g = 2), will have the same possible fates, but with new probabilities (  ( = 2) and   ( = 2)).If the cell enters its PRL, it will have the probabilities to either die in the next generation (  ( = 2) ) or stay in its post-replicative lifespan (  ( = 2)).(C) Probability functions in all conditions.At the beginning of a cell's lifespan, the cell has a high probability to divide, and a low probability to die or enter its post-replicative lifespan.The probability to enter PRL increases during the middle of its lifespan, and the probability of death increases during the end of its lifespan.

Fig. S8 .
Fig. S8.Distribution of lifespans in terms of number of generations and survival functions for different phenotypes in three different media (A) The entire lifespan (aRL + PRL) of each phenotype is represented for two strains (wt and ΔatpA).(B) The apparent replicative lifespan (aRL) of each phenotype is represented for two strains (wt and ΔatpA) in terms of the number of generations.(C) Survival functions of wild-type and ΔatpA strains are divided into two phenotypes.The dotted lines are the experimental data, the solid lines are fits to a Gompertz equation, where (ℎ −1 ) is the aging rate presented on the graph and the numbers in the parentheses are 95% confidence bounds.Red represents phenotype I and blue represents phenotype II.To account for the uncertainty of whether the surviving cells at the end of the experiment are phenotype I or II, we excluded them from the analysis, resulting in all survival functions to go to zero at the experiment's endpoint.

Table S1 .
15st of isogenic single-gene-deletion E. coli K-12 strain BW25113 mutants used in this study.All strains were obtained from the Keio collection.15

Table S2 . Filamentation and death phenotype ratios of individual ∆atpA and wt cells grown in three growth media (LB, M9CG and M9G).
Exp. #. indicates individual experiments.The total number of cells in this table indicates how many cells were analyzed to determine the percentage of the phenotypes.This value might be different of that presented in Figure Movie S1 (separate file).Movies of cell growth, division, and death for example wt cells of phenotypes Ia, Ib, and II.The images are an overlay of DIC, the chromosome in green, and propidium iodide in red.