Spurious regulatory connections dictate the expression-fitness landscape of translation termination factors

During steady-state cell growth, individual enzymatic fluxes can be directly inferred from growth rate by mass conservation, but the inverse problem remains unsolved. Perturbing the flux and expression of a single enzyme could have pleiotropic effects that may or may not dominate the impact on cell fitness. Here we quantitatively dissect the molecular and global responses to varied expression of translation termination factors (peptide release factors, RFs) in bacterium Bacillus subtilis. While endogenous RF expression maximizes proliferation, deviations in expression lead to unexpected distal regulatory responses that dictate fitness reduction. Molecularly, RF depletion causes expression imbalance at specific operons, which activates master regulators and detrimentally overrides the transcriptome. Through these spurious connections, RF abundances are thus entrenched by focal points within the regulatory network, in one case located at a single stop codon. Such regulatory entrenchment suggests that predictive bottom-up models of expression-fitness landscapes will require near-exhaustive characterization of parts. Highlights Precision measurements enable multiscale expression-to-fitness mapping. RF depletion leads to imbalanced translation for co-transcribed gene pairs. Imbalanced translation induces unintended regulons to the detriment of cell fitness. Swapping a single stop codon rewires global susceptibility to RF perturbation.


Introduction
Formulating predictive models connecting genetic information to phenotypes constitutes an overarching goal in genomics and systems biology (Lopatkin and Collins, 2020;Ostrov et al., 2019;Shendure et al., 2019). In single-celled microbes, the relationship between genotype and phenotype can be conceptually decomposed into two distinct maps; the first relating genome sequence to gene expression, and the second, connecting gene expression to whole-cell properties such as proliferation. Rapid progress on the characterization of cis-regulatory elements, spurred by integration of massively parallel reporter assays (Patwardhan et al., 2009(Patwardhan et al., , 2012Sharon et al., 2012) with novel computational frameworks (de Boer et al., 2020;Bogard et al., 2019;Jaganathan et al., 2019;Rosenberg et al., 2015), has achieved headway in predicting expression features from DNA sequence (Agarwal and Shendure, 2018;Cambray et al., 2018;Sample et al., 2019;Urtecho et al., 2020). By contrast, expression-fitness landscapes, defined as the distinct relationships between the expression level of individual genes and the cell fitness, remain understudied despite being the basis of selective pressures on protein abundances. This information gap limits both the engineering of complex biological functions and the interpretability of genetic variation.
Predicting the shape of expression-fitness landscapes requires quantitative characterization of the cellular state at multiple levels. Although numerous expression-fitness landscapes have been previously mapped (Chou et al., 2014;Dekel and Alon, 2005;Duveau et al., 2017;Hawkins et al., 2019;Jost et al., 2020;Kavčič et al., 2019;Keren et al., 2016;Knöppel et al., 2016;Li et al., 2014;Palmer et al., 2018;Parker et al., 2020;Schober et al., 2019;Tubulekas and Hughes, 1993), these measurements rarely include concomitant assessment of the internal cell state following perturbations (but see (Jost et al., 2020;Parker et al., 2020)). With limited information bridging the molecular to cellular scales, the root causes of observed fitness defects are challenging to pinpoint (Fig. 1a). In particular, changes in enzyme levels not only directly affect flux and growth (Fig. 1a, inset i for the case of translation factors) (Ehrenberg and Kurland, 1984;Klumpp et al., 2013;Li et al., 2014), they can also have indirect pleiotropic effects that take the form of damage propagation across three connected levels of biological organization (Fig. 1a, inset ii). First, a reduced enzymatic flux could affect the expression of other genes via specific molecular mechanisms (mechanistic level). Second, these proximal changes in expression could ripple through the regulatory network, leading to further changes in expression genome-wide (regulatory level). Third, each terminal downstream expression change could have an impact on fitness (systemic level). Whether selective pressures on the abundance of proteins predominantly relate to direct impacts on flux or are rather dominated by indirect cellular responses remains unresolved.
Here as a case study, we systematically vary the expression of enzymes involved in the core process of mRNA translation termination (Fig. 1b) in Gram-positive bacterium Bacillus subtilis. We focus on peptide chain release factors, RF1, RF2, and their methyltransferase PrmC (hereafter collectively referred to as release factors, RFs). RF1 and RF2 catalyze the first step of translation termination (Bertram et al., 2001), recognizing stop codons and releasing completed peptides from the ribosome (Fig. 1b). RF1 and RF2 directly interact with the ribosome (Capecchi, 1967;Capecchi and Klein, 1970;Petry et al., 2005), and have partially overlapping specificities (Scolnick et al., 1968) (RF1 recognizes stops UAA/UAG, and RF2 stops UAA/UGA). PrmC post-translationally modifies RF1 and RF2, thereby increasing their catalytic activity (Heurgué-Hamard et al., 2002;Mora et al., 2007;Nakahigashi et al., 2002). We targeted translation termination because the resulting downstream changes in expression were anticipated to be modest: given that translation is initiation-limited on most mRNAs (Laursen and Sørensen, 2005), mildly decreasing termination rate was expected to reduce global ribosome availability without altering protein production on a gene-by-gene basis. A better understanding of the physiology of translation termination stress has relevance in synthetic biology, for example in the context of genome-wide stop codon reassignment (Fredens et al., 2019;Johnson et al., 2011Johnson et al., , 2012Lajoie et al., 2013;Wannier et al., 2018).
Through precise measurements of transcriptomes, global translational responses, and cell fitness (defined here as the population growth rate in exponential phase), we elucidate underlying determinants of the RF expression-fitness landscapes (Fig. 1c). We find that idiosyncratic and indirect inductions of regulatory programs are associated with decreases in growth rate in multiple directions of the RF expression subspace. We term such situation 'regulatory entrenchment,' whereby the fitness defect caused by perturbing a protein's expression is strongly and spuriously aggravated by the gene regulatory network. Further, we reconstruct links that connect the initial microscopic perturbation to system-wide effects. At the mechanistic level, we find that occlusion of ribosome binding sites during RF depletion is a common phenomenon leading to changes in expression stoichiometries between adjacent, cotranscribed genes. In particular, we identify the stop codon of a single regulator as a molecular Achilles heel, sensitizing the entire cell to specific RF perturbations. At the regulatory level, we show that removing one gene can mute pleiotropic changes and liberate RF from regulatory entrenchment. Finally, at the systemic level, we report passive proteome compression as a quantitatively tractable cause of growth defects upon massive activation of a regulon.
Our multiscale characterization provides a concrete example of how distal yet focal events triggered by targeted molecular perturbations can have system-wide impacts. Here, the cells' susceptibility to expression perturbations of specific enzymes is not simply related to the magnitude of the changes in the associated cognate flux, but is instead dictated and amplified by sensitive nodes in the regulatory network. These results underscore the viewpoint that a quantitative understanding of the full system (Boyle et al., 2017;Karr et al., 2012;Liu et al., 2019) might be necessary to predict even qualitatively the shape of expression-fitness landscapes.

Linking RF perturbations to changes in genome-wide expression and fitness
We used a combination of genetic tools and high-throughput measurements to map RF expression-fitness landscapes, as well as the underlying mechanistic, regulatory, and systemic responses. We created strains in which two of the three factors (PrmC, RF1, and RF2) can be tuned orthogonally (Fig. S1e-f): PrmC and RF1 in one strain, and PrmC and RF2 in another (with the autoregulatory frameshift removed for RF2 (Craigen and Caskey, 1986;Craigen et al., 1985), Fig. S1a-c). The range of these tunable expression systems spanned 31 to 111-fold and was centered near their respective endogenous levels (Fig. S1d). The global gene expression changes resulting from these perturbations were probed with RNA-seq (Parker et al., 2019) (Methods). Expression levels were converted to fractions of the proteome using ribosome profiling (Ingolia et al., 2009;Li et al., 2014) (Methods for details and assumptions). Ribosome profiling further provided a high-resolution view of translation status in a subset of conditions. In order to precisely measure fitness, defined here as relative population growth rate in exponential phase, we designed a DNA-barcoded competition assay (Parker et al., 2020;Smith et al., 2009) with ±1% precision (±2σ, Fig. S2, Methods). This combination of targeted proteome perturbation with precision transcriptomic and fitness measurements enabled a multiscale assessment of the RF expression-fitness landscapes.

Perturbing the expression of different RFs decreases fitness through distinct physiological routes
Using our measurement platform, we directly confirmed that endogenous RF expression in exponential phase maximizes growth rate. For conditions at or near endogenous levels of RF1, RF2, and PrmC (dashed lines, Fig. 2e-g), the fitness is indistinguishable from wild-type strains (shaded grey area in Fig. 2e-h corresponds to experimental precision in fitness, |s| < 2σs = 1.2%, Methods). On the other hand, perturbing RF expression away from endogenous level led to growth defects in most directions of the RF subspace ( Fig. 2e-h). Only RF1 overexpression did not cause a measurable decrease in fitness, in part because of the limited maximal achievable level for this particular inducible construct (≈3×, Fig. S1d). These results suggest that RF expression is optimized for exponential phase growth, akin to the expression of several other factors (Dekel and Alon, 2005;Ehrenberg and Kurland, 1984;Li et al., 2014;Parker et al., 2020), and is consistent with the tight evolutionary conservation of expression stoichiometry in biological pathways (Lalanne et al., 2018).
Away from the optimum, distinct RFs displayed expression-fitness landscapes of different shapes. For example, PrmC has a much more severe growth defect than RF2 at the same levels of overexpression ( Fig. 2f vs. 2g), and RF1 knockdown leads to a near-vertical drop in fitness (Fig. 2e), which is more severe than RF2 or PrmC knockdown (though our expression constructs had higher baseline expression for the latter two). Beyond these qualitative differences, our global expression quantification in all these conditions provides a way to assess whether RF-specific fitness defects nevertheless correspond to similar underlying physiological states generic to translation termination defects.
Projecting complex phenotypic data onto low-dimensional manifolds can provide insight on the physiological state underlying growth defects (Hui et al., 2015;Scott et al., 2010Scott et al., , 2014You et al., 2013). Two cellular-level quantities are of particular interest: (1) the growth rate , and (2) the total expression of the mRNA translation machinery, quantified as the summed proteome fraction of translation proteins, termed the translation sector, R. Hwa and coworkers have shown (Scott et al., 2010;You et al., 2013;Zhu et al., 2016Zhu et al., , 2019) that trajectories in the space of growth vs. translation sector R upon a series of increasingly severe perturbations are intimately related to global regulation (Box. 1). In E. coli, they found that under numerous ways to inhibit translation, the translation sector R increases concomitantly with a decrease in the growth rate (along the "translation line", Box 1) as a compensation mechanism (Scott et al., 2010). By contrast, decreasing the nutritional quality of the medium leads to a decrease in growth paralleled by a decrease in the translation sector R (along the "nutrition line", Box 1), a long-known property in bacterial physiology (Schaechter et al., 1958). Under the bacterial growth laws established in E. coli (Box. 1), the states of cells with RF-specific fitness defects were anticipated to collapse on the translation line. In effect, perturbing RF expression was expected to reduce the translation termination rate and thus to reduce ribosome availability.
Surprisingly, the different RF expression perturbations displayed diverse physiological trajectories ( Fig. 2i-l). The sharp drop in growth under RF1 knockdown led to little change in the translation sector R (Fig. 2i), analogously to the trajectories associated with non-translation targeting antibiotics (Scott et al., 2010). By contrast, the growth defects upon PrmC expression perturbations were associated with a decrease in the translation sector ( Fig. 2k-l), leading to movement along the nutrition line, which was unexpected given the fixed quality of the growth medium. RF2-perturbed cells moved slightly in different directions following overexpression (translation line) and knockdown (nutrition line, although high basal activity of our expression construct limited the magnitude accessible growth defects) (Fig. 2j). Interestingly, physiological changes were additive under combined RF perturbations: tuning PrmC expression in combination with RF2 overexpression led to movement along the nutrition line, but shifted by the movement along the translation line by the RF2 perturbation (dashed vs. dotted lines in Fig. 2j-l, Supplementary Discussion).
The physiological divergences observed for the different RFs perturbations were intriguing given the similar involvement of these three proteins in translation termination, and pointed to possible RF-specific pleiotropic effects. Analyzing the transcriptomes following perturbations revealed distinct responses along the orthogonal RF expression directions. RF1 knockdown led to changes in motility and biofilm genes ( Fig. S4c-d), PrmC perturbations massively induced the general stress B regulon ( Fig. 3a-b, S3a-d), and RF2 perturbation led to little changes across the full range of expression, except for a modest B induction at maximal knockdown ( Fig. S3o-p). Together, these results suggest that each RF has a mechanistically distinct relationship between expression and fitness.

Fitness defects during PrmC overexpression can be traced back to regulatory response by B activation
To further identify the origin of some of the observed growth defects, we focused on perturbations to PrmC. The growth defect following PrmC overexpression ( Fig. 2g) coincided with dramatic activation of the general stress regulon, which is driven by the alternative sigma factor B (gene sigB) (Hecker et al., 2007;Helmann et al., 2001;Petersohn et al., 2001;Price et al., 2001;Price CW, 2002;Zhu and Stülke, 2018), from 2% basal proteome synthesis fraction in wild-type to >20% (Fig. 3a, c, see Methods for ribosome profiling calibration). This increase in the expression of B genes was accompanied by a corresponding decrease in the translation sector R (Fig. 3d-e), which raised the possibility that a substantial portion of the growth defect following PrmC overexpression was associated with the physiological burden of expressing >100 B -dependent genes at high levels (Benson and Haldenwang, 1992;Bernhardt et al., 1997;Boylan et al., 1992). Alternatively, induction of the stress regulon might have directly responded to, and played a role in mitigating, translational stress resulting from PrmC expression changes.
As additional evidence of systemic stress, PrmC overexpression leads to global changes in gene expression beyond the B regulon (broadening of fold-change distribution, Fig. 3a inset). Whether these global changes were the result of the PrmC perturbation, or secondary to B activation, was to be clarified.
To determine whether B exacerbates or mitigates fitness defects under PrmC perturbations, we profiled RF-inducible cells lacking the sigB gene (∆sigB). As expected, deleting sigB completely abrogated regulon induction ( Fig. 3b-c, dashed red arrows), but surprisingly also restored genome-wide expression to near endogenous levels ( Fig. 3b inset). Hence, global expression changes (gray line in inset Fig. 3a) upon PrmC overexpression were caused by B activation and its downstream consequences (e.g., sigma factors competition (Farewell et al., 1998)). Importantly, the PrmC overexpression growth defect was completely rescued following sigB deletion ( Fig. 3g dashed red arrow), concomitantly with the restoration of the translation sector to endogenous levels ( Taken together, these results suggest that the gene sigB constitutes a trans-modifier (Hou et al., 2019) of the PrmC expression-fitness landscape. Through its activation, B drastically exacerbates the cell's sensitivity to expression perturbation of distal gene prmC. Overall, the induction of the general stress regulon has no measured benefit in the context of PrmC perturbations. Gene sigB therefore 'entrenches' PrmC expression by imposing a massive physiological burden to cells that is in addition to, and much larger than, the effect directly related to perturbing PrmC.

Systemic proteome compression by B explains fitness defect
We hypothesized that the growth defect resulting from B activation mainly arose from the cost of expression of regulon genes. The burden of expressing proteins has been extensively characterized in E. coli (Andrews and Hegeman, 1976;Dong et al., 1995;Scott et al., 2010;Stoebel et al., 2008). Gratuitous expression of genes causes "proteome compression", whereby the abundances of other proteins globally decrease in response to the finite total synthesis capacity of the cell (Scott et al., 2010). The resulting reduction in the translation sector, responsible for generating biomass, leads to a growth defect, which would be in addition to any protein-specific effects such as misfolding or aggregation. Bacterial growth laws, established in E. coli, predict the quantitative relationship between proteome fraction of gratuitous proteins and growth rate under this proteome compression model (Box 1).
Despite the complexity of the induced regulon (>100 diverse genes) and the resulting secondary changes in gene expression (Fig. 3a), B activation did lead to simple compression of the translation machinery ( Fig. 3e-f). Indeed, the translation sector (full line, Fig. 3e, S3g) and growth rate (full line, Fig. 3h, S3i) linearly decreased in a one-to-one proportion with the increase in excess proteome fraction to the B regulon. Interestingly, we observed a smaller growth defect than expected from the E. coli growth law (dashed vs. full lines, Fig. 3e, 3h, S3g, S3i, Box 1). Further assessment of the generality of growth laws to species other than E. coli will constitute promising future research avenues.
Compression of the translation sector due to B activation provides a systemic explanation for the movement along the nutritional line under PrmC perturbation (Fig. 2k-l). Abrogation of this effect following sigB deletion ( Fig. 3d-e, S3f-g, dashed arrow) highlights that physiological trajectories can be modulated by individual regulatory proteins, as has been also demonstrated for tRNA synthetases and the stringent response in B. subtilis (Parker et al., 2020). However, unlike the stringent response, which is cognate to synthetase perturbation and beneficial to cell growth, the B response here is detrimental.

Activation of regulator B does not result from general stress-sensing mechanisms
The detrimental effect on fitness of B activation under RF stress (Fig. 3g) suggests that its induction is idiosyncratic in this context, as opposed to being the result of a cognate sensing mechanism. We hypothesized that RF-specific perturbations led to changes in of a small number of proteins that ultimately activated B . In support of this view, B is only activated when certain RFs are perturbed: varying RF1 expression over the full achievable range (0.1× to 3× endogenous) did not change B regulon expression (Fig. S4c), whereas RF2 knockdown, but not overexpression, was associated with B activation (Fig. S3o-p, and lowering basal RF2 expression further increased B activation, see Supplementary Data 8). Similar to PrmC overexpression, the fitness defects under RF2 knockdown was rescued by sigB deletion (Fig.  S3q), also arguing for regulatory entrenchment of RF2 expression by B .
To confirm that B activation arose from perturbations specific to distinct RFs and not from general translation stress-sensing mechanisms, we used ribosome profiling to quantify translation in cells acutely depleted of RFs. We used CRISPRi transcriptional interference (Peters et al., 2016;Qi et al., 2013) to separately knockdown RF1/PrmC (co-transcribed) and RF2. In both cases, we found evidence of translational stress in the forms of idle ribosomes at the corresponding stop codons and queuing upstream (Fig. 4a), similar to what was previously observed under different translation termination stresses (Baggett et al., 2017;Mangano et al., 2020;Saito et al., 2020). Specifically, we observed queues at UAG stops under RF1/PrmC depletion, queues at UGA stops under RF2 depletion, and no queues in wild-type, or at UAA stops in either depletion. Queues were longer on mRNAs with high translation efficiency, as expected from models of stochastic queuing (Baggett et al., 2017;Bergmann and Lodish, 1979;Mitarai et al., 2008;Saito et al., 2020) (Fig. S5b-c, S5h, Supplementary Discussion). The severity of the translation stress, as assessed by the magnitude of accumulation of ribosomes on stop codons, was similar for the PrmC/RF1 and RF2 depletions (Fig. 4a). Importantly however, a robust B activation was only observed under RF2, and not RF1/PrmC, knockdown (Fig. 4b-c). The simultaneous presence of ribosome queues (Fig. 4a, middle row) and absence of B activation under RF1/PrmC depletion (Fig. 4b) confirmed that B induction was not channeled through a sensor of translation stress agnostic to the identity of ribosome-jammed mRNAs. B induction was instead likely due to gene-specific expression changes driven by RF stop codon specificities.

RF depletion mechanistically leads to stoichiometric imbalance of co-transcribed genes
We hypothesized that a consequence of RF depletion is that the ribosomes idling at the corresponding stop codons would sterically occlude translation initiation of downstream genes. Because bacterial genes are often closely spaced in polycistronic operons, elevated ribosome occupancy at stop codons may prevent other ribosomes from binding to the next gene (Fig. 4d). Consistent with this hypothesis, we found that among co-directional genes separated by a ribosome footprint or less (≤30 nt), the downstream genes exhibit a mild but stop-specific decrease in expression when the cognate RF is depleted (Fig. 4d, RF1/PrmC: FCUAG = 0.88, p<10 -6 ; RF2: FCUGA = 0.82, p<10 -6 , see Fig. S5d-f for control groups). This result differs from a lack of effects reported in E. coli during knockdown of ribosome recycling factors (RF4), which act on ribosomes post RF1/RF2-mediated peptide release. These ribosomes may have different mobility and re-initiation properties (Saito et al., 2020). Although the effect size was modest overall, some gene pairs were affected more heavily. For example, highly translated upstream genes are correlated with stronger decrease in downstream expression Fig. S5i. Gene pairs joined by AUGA, i.e., the most common overlapping start and stop punctuation (Fig. S5j) possibly involved in translational coupling, did not show a substantial difference (Fig. S5k). Overall, these genome-wide measurements point to occlusion of ribosome binding sites by upstream termination-idle ribosomes as a cause of RF-specific gene expression changes.
In contrast to the effect on co-transcribed genes downstream, the expression of genes ending with the compromised stop codons were not substantially affected either at the level of translation (Fig. 4e . This lack of systematic changes in expression for genes with translation termination defects is consistent with the observed ribosome queues being too short to interfere with translation initiation at the upstream start of the ORF (queue size of <4 stalled ribosomes globally, Fig. 4a) (Bergmann and Lodish, 1979;Jin et al., 2002). This result also suggests that translational quality control and surveillance mechanisms, such as mRNA cleavage at empty A-site (Hayes and Sauer, 2003;Li et al., 2007) and ribosome collision sensors (Ferrin and Subramaniam, 2017), play a minor role in comparison to start codon occlusion for downstream genes in the current context.

Mechanism for B activation by hypersensitivity of a single cis-element
Expression stoichiometry of genes rsbW/rsbV, whose ORFs overlap via the RF2 stop codon (-4 nt overlap AUGA, Fig. 5b, Fig. S5j), was one of the most sensitive to RF2 knockdown (3.0× change, Fig. 4d, compared to 1.1× in RF1/PrmC knockdown). Intriguingly, these two proteins form the regulatory core of B activation. Upstream signaling events converge on a phosphorylation-dependent partner-switching mechanism involving the anti-sigma factor RsbW and the anti-anti-sigma factor RsbV, which respectively inhibits and activates B (Pané-farré et al., 2017). As additional control points, RsbW can deactivate RsbV by direct phosphorylation (Hecker et al., 2007;Price CW, 2002) and the three genes rsbW, rsbV, and sigB are cotranscribed under a B dependent promoter (among other transcript isoforms) (Fig. 5a). These interlocked positive and negative feedbacks (Fig. 5a) render B activation hypersensitive to the stoichiometry of its regulators RsbV and RsbW (Igoshin et al., 2007;Locke et al., 2011).
Given that we observed B activation in RF2, but not RF1/PrmC, knockdown (Fig. 4c  vs 4b), and that the expression of B regulators RsbV/RsbW was affected in these conditions (Fig. 4d), we hypothesized that termination defects at the rsbV UGA (RF2-specific) stop codon was the cause of the B activation following RF perturbations. As such, switching the rsbV stop codon was predicted to be sufficient to rewire the susceptibility of B activation to different RF perturbation. For example, with a switched UAG stop for gene rsbV, B might be activated upon RF1, but not RF2, knockdown.
To test whether B activation sensitivity was indeed focally dependent on this ciselement, we generated a set of rsbV stop codon variants (UAA: RF1 and RF2 activity, UAG: RF1 only activity) (Fig. 5b). These variants were cleanly inserted at the endogenous sigB operon in wild-type and in an RF inducible strain (RF1/PrmC and RF2, with frameshift removed, under orthogonal promoters, Methods). Under various RF expression conditions (endogenous, RF2 knockdown, RF1/PrmC knockdown), we monitored mRNA levels of two responsive Bdependent genes by RT-qPCR (ywzA and ygxB, highlighted in Fig. 3a, Methods). Basal B activity at endogenous RF expression differed slightly for the different rsbV stop variants, likely as a result of perturbed ORF positioning. Strikingly, we observed strong B activation only in strains with the rsbV stop cognate to the RF perturbation (Fig. 5c, raw data in Supplementary Data 8). Consistently with previous experiments, RF2 knockdown led to an increase (>8×) in the expression of our B reporter genes in the rsbVUGA variant (Fig. 5c, top row). In addition, induction of reporter genes was observed for the rsbVUAG variant only under RF1/PrmC knockdown (>50×, Fig. 5c, bottom row). No induction in either RF perturbation was observed for the rsbVUAA variant (Fig. 5c, middle row), in accordance with the UAA stop codon being RF agnostic. The rsbVUAA variant was confirmed functional (>500× induction, indistinguishable from the wild-type rsbVUGA variant following 5 minutes of 4% v/v ethanol stress, Supplementary Data 8). Hence, our stop codon swapping experiment demonstrates that systems-wide susceptibility to a distal molecular perturbation can be encoded through a single sensitive ciselement, and that such susceptibility can be entirely rewired with minimal genetic changes.

Discussion
We have shown that quantitative genome-wide measurements in conjunction with precision fitness quantification can be used to reconstruct chains of events across biological scales relating a molecular perturbation (change in expression of a specific protein) to a whole cell phenotype (growth rate). We find that in B. subtilis, important portions of the selective pressures on peptide release factor abundances is through idiosyncratic activation of distal endogenous regulatory programs (Fig. 6). By analogy with interactions between residues constraining the evolution of individual proteins (Bridgham et al., 2009;Pollock et al., 2012;Shah et al., 2015;Starr et al., 2018), we term this strong dependence of enzymes' expression-fitness landscapes on the network of genetic interactions 'regulatory entrenchment' (Fig. 6).
Despite our focus on translation termination factors, which we anticipated to have limited impacts on gene expression (translation termination being non-limiting on mRNAs), pleiotropic effects related to activation of specific regulons still dominated the connection between RF expression and fitness. As an example of molecular determinants shaping expression-fitness landscapes, we found that a single cis-element (identity of the rsbV stop codon, Fig. 5b-c) mechanistically led to RF-specific expression imbalance (altered stoichiometry of B regulators RsbW and RsbV, Fig. 4d). This imbalance was amplified to global transcriptome remodeling through to the activation of B (Fig. 3a, S3c). These changes ultimately had systemic impacts (growth defects by proteome compression, Fig. 3e, S3g, S3o). Although the exact cause of B activation due to PrmC overexpression remains to be determined, deletion of the sigB gene was sufficient to flatten the landscape in two directions of the RF expression space (PrmC overexpression Fig. 3g, S3h and RF2 knockdown Fig. S3q).
The adaptive nature of B activation upon RF-stress remains a possibility, but the lack of corresponding fitness advantage in our experiments suggests otherwise (Fig. 3g, S3h). Bioinformatic analysis of homologous B operons shows some conservation in the rsbV/rsbW ORF overlap (Fig. S6, Methods). Overlap of these genes could however be driven by the requirement of a precisely tuned expression stoichiometry (Igoshin et al., 2007;Locke et al., 2011), enabled by translational coupling via the -4 nt start/stop overlap configuration AUGA. In such scenario, susceptibility to RF perturbation would emerge as an evolutionary byproduct of selective pressures that operate on other features of the system (Gould, 1978).
As an additional example of regulatory entrenchment in our data, the sharp decrease in fitness upon RF1 knockdown (Fig. 2e) was associated with gene expression changes characteristic of the motile-to-sessile bistable switch in B. subtilis (Chai et al., 2009;Kearns and Losick, 2005) (Fig. S4). Given the complexity of the regulation of motility genes (Mukherjee and Kearns, 2014), we have not been able to identify a unique plausible molecular link to RF1 knockdown (Supplementary Discussion). The relationship between the transition to sessile state and the plummeting cell growth rate also remains elusive. These underscore the challenges of reconstructing causal features underlying expression-fitness landscapes, even with detailed under-the-hood information about the cell's internal state. Our case study highlights that exhaustive information at multiple levels (Fig. 1a, inset ii), down to the response of all cisregulatory elements, might be required to achieve fully predictive models given how easily perturbations propagate across scales in biological systems. A lack of global characterization of the susceptibilities and interaction of components in regulatory networks currently constitute an important bottleneck in systems biology, although at-scale methods are being developed to dissect genetic networks (Calderon et al., 2020;Muller et al., 2020).
We anticipate that master regulators will be common focal points of regulatory entrenchment, as seen here with general stress factor B and motility factor D , due to the responsiveness of their activation and their direct impact on a large number of genes. As such, successful predictions will presumably be more easily achieved in near-minimal organisms with limited regulation (Baby et al., 2018;Hutchison et al., 2016). Given the broadly conserved expression stoichiometry of enzymes across evolution (Lalanne et al., 2018), we speculate that endogenous expression levels satisfy generic optimization principles, such as flux maximization under resource allocation constraints (Kurland and Ehrenberg, 1987). However, since diverse bacteria harbor distinct cis-elements, operonic organization, and regulatory architectures, the shape of expression-fitness landscapes away from the optima might generally differ as a result of regulatory entrenchment.

Fig. 1. Mapping the underlying determinants of the release factor expression-fitness landscapes. (a)
Expression-fitness landscapes, which connect enzyme expression (microscopic variable) to the growth rate (cellular phenotype), can be dictated by direct, or indirect effects. Inset (i): direct effects correspond to reduction in the flux cognate to the perturbed enzyme (protein synthesis rate in the case of translation factors). Inset (ii): indirect effects result from pleiotropic propagation across mechanistic, regulatory, and systemic levels. (b) As a case study, the expression of peptide chain release factors (RFs: RF1, RF2, and associated methyltransferase PrmC), involved in the first step of mRNA translation termination, was tuned around endogenous levels. (c) Strains with inducible copies of RFs, and deleted endogenous genes, were used to systematically vary RF expression. The resulting impacts on the cell internal state (RNA-seq, ribosome profiling) and relative growth rate s (competition experiments) were measured, leading to precise mapping of expression-fitness landscapes. See Fig. S1 and S2 for details on strains and measurement platform.      Fig. 6. Regulatory entrenchment in gene expression-fitness landscapes. Perturbation of specific enzyme expression can lead to changes in the expression of a subset of genes through the sensitivity of cis-regulatory elements. These expression imbalances can further reverberate in the network of regulatory interactions in the cell through focal points such as master regulators, which control the expression of multiple genes (purple box). In the absence of these regulators, downstream responses are muted (teal box). Pleiotropic expression changes exacerbate the fitness defect of perturbing the original enzyme. The cell's susceptibility to expression is then not simple relation of the decreased enzymatic flux, but also depends on the aggregated impacts regulatory susceptibilities of the whole cell. Cells can be liberated from regulatory entrenchment by removing nodes in the network. See also Fig. S3

Box 1. Bacterial growth laws
The bacterial growth laws (Scott et al., 2010(Scott et al., , 2014You et al., 2013) constitutes an empirical framework recapitulating the global regulatory architecture of E. coli under nutrition and translation stress in steady-state growth. The growth laws connect the growth rate to the abundance of coarse-grained proteome sectors (translation or "ribosome affiliated" proteins + ; unregulated proteins which includes metabolic enzymes 2 , and gratuitous unnecessary non-toxic proteins 3 ) through three mathematical relationships: Resource allocation: + ,-.
Above, phenomenological parameters ∘ corresponds to an inactive ribosome fraction, ' is the nutritional capacity, ) is the translational capacity (proportional to the rate of protein synthesis), and is a fixed conversion factor. Nutritional and translational capacities are defined through slopes of the lines along physiological trajectories ( Fig.  Box. 1).
The growth laws can be interpreted as amino acid flux conservation in steady-state growth: the flux generated by catabolic enzymes and transporters ( 2 ) matches the biomass production flux by translation proteins ( + ). Changing the medium quality changes nutritional capacity ' without altering the translational ) and inactive ribosome content ∘ . Then, the translation sector + decreases following the nutrition line: Nutrition line: + = ∘ + ) .
Upon expression of unnecessary proteins to proteome fraction 3 , the growth laws above predict movement along the nutrition line (since ) and ∘ are assumed fixed), and a fractional decrease in the translation sector and growth: The slope of + vs. 3 is determined by, which relates to the posited incompressibility of fraction 1 − + ,-. of the proteome. The different slopes (Fig. 3e, 3h, S3g, S3i) we observe in B. subtilis suggest possible differences in the growth laws in this species.