Multiscale effects of heating and cooling on genes and gene networks.

Most organisms must cope with temperature changes. This involves genes and gene networks both as subjects and agents of cellular protection, creating difficulties in understanding. Here, we study how heating and cooling affect expression of single genes and synthetic gene circuits in Saccharomyces cerevisiae We discovered that nonoptimal temperatures induce a cell fate choice between stress resistance and growth arrest. This creates dramatic gene expression bimodality in isogenic cell populations, as arrest abolishes gene expression. Multiscale models incorporating population dynamics, temperature-dependent growth rates, and Arrhenius scaling of reaction rates captured the effects of cooling, but not those of heating in resistant cells. Molecular-dynamics simulations revealed how heating alters the conformational dynamics of the TetR repressor, fully explaining the experimental observations. Overall, nonoptimal temperatures induce a cell fate decision and corrupt gene and gene network function in computationally predictable ways, which may aid future applications of engineered microbes in nonstandard temperatures.

Temperature plays a critical role across all of biology. Endothermic organisms, commonly called “warm-blooded,” maintain their body temperature and perish if this regulation fails. Despite common portrayal of most species as “cold-blooded,” even ectotherms and unicellular microbes prefer a specific temperature range where they function optimally. If the ambient temperature departs from this optimum, cells and organisms attempt to minimize the potential damage. For instance, the heat shock response elicits the transient expression of cytoprotective proteins, mitigating a multitude of harmful effects ranging from changes in cell morphology to alteration in protein structure and function (1–3). Lack or failure of heat shock response leads to cell cycle arrest, or sometimes cell death (4–6).

Gene-regulatory networks govern cellular protein levels and play a double role during temperature changes. First, they require protection due to their crucial role in cellular housekeeping and homeostasis (7), development (8), and survival (9). Second, they are also instrumental in generating the protective stress response. This twofold involvement as both the subject and agent of cellular protection makes it difficult to understand what happens with gene networks at nonoptimal temperatures. Additional complexity stems from nontrivial effects of temperature on transcription (10), translation (11), and their regulation (12, 13). Previous transcriptome expression measurements in heat and cold shock revealed broad responses (14), but it is hard to discern the degree to which these changes are individual gene-level (15) or network-level effects, or cellular attempts to restore homeostasis. Adding to these difficulties are the complexity of biochemical interaction networks and the incomplete knowledge of their connectivity (16, 17).

Synthetic gene circuits are relatively simple human-designed gene networks that perform specific predefined functions, and have promising roles as switches (18), oscillators (19, 20), and logic gates (21) in future medical (22–24), bioenergetic (25), or environmental (26) applications. The connectivity and components of synthetic gene circuits are typically well characterized and completely known by design, and their regulatory connections with the host network are assumed minimal. Therefore, these relatively simple gene network modules should be less affected by cellular attempts to restore homeostasis, making them viable candidates for discerning and separately studying the network module-level effects of temperature changes. Besides a potential to advance fundamental understanding of temperature effects on biological systems, future synthetic biology applications will require accuracy and robustness in carrying out predefined functions, independent of environmental changes (27, 28). In some cases, temperature provides a way to control synthetic gene circuits (29) and temperature compensation may be possible (30). However, since most synthetic gene circuits are developed and characterized under controlled laboratory conditions, it is unclear how temperature shifts affect their function (31), such as the gene expression noise they are designed to impose.

To unravel how temperature affects genes and gene networks, we chose to study multiple single genes and two synthetic gene circuits: the “negative-feedback” (NF) or “linearizer” (32) and “positive-feedback” (PF) (33) gene circuits. Positive and negative autoregulation are highly abundant motifs in natural gene networks (34–36), and inducible systems with autoregulation are common in synthetic biology (37–39). Since the PF and NF gene circuit components are not native to yeast and lack substantial connections with the native yeast regulatory network, their behavior at 30 °C was predictable by quantitative modeling (32, 33). Therefore, any deviations from this standard behavior are potential nonoptimal-temperature effects. To identify such effects, we measured both PF and NF gene expression and cell division rates, which mutually influence each other (33, 40–42). Apart from suboptimal growth below and above 30 °C, we discovered unusual gene expression bimodality originating from a critical cell fate decision between stress resistance and growth arrest with vanishing gene expression. The resistant cells had nontrivial gene expression trends. Computational models with temperature-dependent growth rates and Arrhenius reaction kinetics captured the functional effects of cooling, but not those of heating in resistant cells. On the other hand, molecular-dynamics (MD) simulations revealed an ability of inducer-bound, normally inactive TetR to repress its target promoter. Once we augmented the quantitative models with this altered protein behavior, we could fully explain the changes in gene circuit function at all temperatures. Our approach of systematically dissecting temperature effects on gene network modules combined with multiscale mathematical modeling and simulations is broadly applicable to analyze temperature effects on other genes and regulatory network modules. So, studying how heating and cooling affect the PF and NF gene circuits in relative isolation from temperature compensation and other natural network complexities gives clues on how heating and cooling generally impact gene network module function. Overall, the cell fate bifurcation, combined with altered growth rates, reaction rates, and protein structure represent four pillars of understanding and predicting the behavior of genetic systems at nonoptimal temperatures in general.

Results

Temperature-Dependent Growth and Expression in Yeast Cells with a Single Reporter Gene.

Since genes and their promoters are the basic building blocks of gene-regulatory networks, we began by investigating how temperature affects a single reporter gene in the YPH500 budding yeast strain (MATα ura3-52 lys2-801_amber ade2-101_ochre trp1-Δ63 his3-Δ200 leu2-Δ1; congenic with the standard laboratory strain S288c). First, we focused on NF0 cells that carried only the reporter gene from the full NF gene circuit (32), but without a regulator ( and Fig. 1). In the “NF0” naming convention, the zero indicates a single gene “without regulation” compared with the strain carrying the complete NF gene circuit. Specifically, the NF0 strain incorporates the pGAL1-based, TetO2 operator site-containing P promoter (43, 44) driving the yEGFP::zeoR reporter gene (Fig. 1), but lacks the tetR repressor.

Fig. 1.

The effect of temperature on NF0 (single reporter strain) growth rates and gene expression. (A) Schematics of NF0 consisting of the yEGFP::zeoR reporter controlled by the P promoter. (B) Full population and resistant-cell (R cell) subpopulation growth rates of NF0 cells as a function of temperature. (C) Representative NF0 gene expression histograms from experiments. (Inset) (100×) An NF0 A cell and an R cell at 38 °C. (D) Fraction of growth-arrested low-expressing NF0 A cells from NF0 gene expression histograms in C. (E) Average reporter (yEGFP::zeoR) expression as a function of temperature from experiments (dark gray) and gene expression results from Arrhenius models (light gray) for NF0 cells (see , for parameters). Temperature-dependent differences are statistically significant (P = 1.586 × 10−13, one-way ANOVA). Gene expression and growth rate values were normalized by the corresponding values from replicates in the 30 °C control condition. Error bars represent the SEM (n = 3). Expt, experiment; norm, normalized.

Considering that cellular growth rates and gene expression mutually affect each other (33, 41, 42), we started out by measuring the growth rates of NF0 cells at temperatures ranging from 10 to 40 °C. The growth rates of NF0 cells had an optimum at 30 °C (Fig. 1 and , which is possibly due to decades-long growth in standard laboratory conditions. Below 30 °C, the growth rates had Arrhenius temperature dependence (R2 > 0.98) (). Above the optimum, growth rates declined rapidly, approximately linearly with temperature.

Next, we studied temperature dependence of gene expression for NF0 cells by flow cytometry. Surprisingly, NF0 gene expression was bimodal at all temperatures above and below 30 °C, in sharp contrast with its normal, unimodal appearance at the optimal temperature of 30 °C (Fig. 1). Past studies have shown that yeast cells can arrest growth at nonoptimal temperatures (45–47) and acquire autophagic signatures (48). Growing NF0 cells at 38 °C for 2 d and then imaging them in a microfluidic chip (49) revealed that some, but not all cells had expanding, granular vacuoles as their nuclei shrank and their fluorescence diminished. This process continued until the vacuole filled such cells, while their fluorescence and then their nucleus completely vanished. These features suggest that some cells might attempt surviving by autophagy (50, 51), until they succumb to cell death. Taken together with the stable bimodal gene expression histograms for cooled and heated NF0 cells, we concluded that cooled and heated yeast cells segregate into two phenotypically distinct subpopulations (Fig. 1, Inset). In arrested, autophagous “A cells,” nonoptimal temperature abolishes growth jointly with gene expression, while resistant “R cells” continue growing and expressing the reporter, maintaining overall growth and stable bimodal distributions at the cell population level.

To examine the origins of this cell fate bifurcation, we reasoned that temperature stress acts as a selective force on a phenotypically heterogeneous, intrinsic cellular property such as Tsl1p (6) or Hsp12p (52) levels. If continued cell growth and gene expression require this property to exceed some threshold, then two subpopulations will naturally emerge (Fig. 2). Specifically, R cells that exceed the threshold will sustain growth and reporter expression, while A cells below the threshold will arrest growth and lose gene expression. Additionally, this simple threshold model predicts that growth in higher stress requires exceeding a higher threshold, which fewer cells can achieve. Consequently, higher stress levels should augment the fraction of A cells at the expense of R cells. Indeed, gene expression measurements confirmed this expectation: The fraction of A cells increased as the temperature deviated from the optimum in either direction (Fig. 1), partly explaining the concave temperature dependence profile of gene expression (). Accordingly, the coefficient of variation (CV) (the ratio of the SD and the mean) for the full NF0 population had a convex trend, with minimum values at 20 and 30 °C ().

Fig. 2.

Temperature effects on gene expression in various yeast strains. (A) Schematic of temperature effects on gene expression and division at the single-cell level, where lighter-gray–shaded R cells are stress resistant, and dark-shaded A cells are sensitive, future nonexpressing, growth-arrested A cells, which eventually become irregularly shaped dark dead cells. At normal temperature (tails of the dashed arrows), cells express proteins and divide. At nonoptimal temperatures (tips of the dashed arrows), only R cells with high levels of resistance (rightward-curving arrows, pointing to R cells with light shading on the Left and green color on the Right) are able to express and divide. This gives rise to the temperature-dependent fraction of growth-arrested A cells (leftward-curving arrows) at the expense of resistant R cells. (B) TBR1 yeast cells have unimodal reporter expression from P at 30 °C, which becomes bimodal at 38 °C. (C) KV2695 yeast cells have unimodal reporter expression from P at 30 °C, which becomes bimodal at 38 °C. (D) Tsl1p-GFP S288c yeast cells have unimodal reporter expression from P at both 30 and 38 °C, as expected for a stress-resistance marker protein. fluor, fluorescence; norm, normalized.

Discounting the nonexpressing A cells, temperature still affected gene expression in the R cell subpopulation (Fig. 1). To examine whether the temperature-dependent trends of gene expression in R cells were purely due to slower growth at nonoptimal temperatures, we developed computational “growth rate models.” These models assumed constitutive gene expression, setting all protein dilution rates equal to R cell growth rates estimated from the experimental temperature-dependent growth rates of the whole population (Fig. 1 and ). These growth rate models could not reproduce the experimentally observed trends for the NF0 R cells (). Therefore, next we also assumed Arrhenius temperature dependence of protein synthesis rates (). These augmented “Arrhenius models” matched the NF0 R cell data well, except at 35 °C (Fig. 1). This discrepancy at 35 °C could arise because the P is a P-based, possibly nonconstitutive promoter, with sugar- and potentially temperature-dependent native regulation. We concluded that Arrhenius models with temperature-dependent growth and reaction rates were sufficient to capture gene expression changes for R cells expressing a single reporter gene at various temperatures. Multiscale models additionally including the population-dynamic effect of growth arrest could explain the full NF0 cell population’s behavior.

Next, we asked whether the temperature-dependent cell fate bifurcation and bimodal gene expression might be observable for some other yeast strains and promoters. Therefore, we investigated how heating affects reporter gene expression from the same P promoter (53) in Σ1278b (also known as TBR1; MATα, ura3-52, his3::hisG, leu2::hisG), which has 3.2 single-nucleotide polymorphisms per kilobase compared with the standard laboratory strain S288c (54, 55). In the TBR1 yeast strain yEGFP gene expression from the P promoter was bimodal at 38 °C and unimodal at 30 °C (Fig. 2), suggesting that the heat-induced cell fate bifurcation is not strain dependent. We also investigated the expression of multicistronic FLO11-2A-YFP (composed of the FLO11 and the YFP reporter genes linked via a 2A self-cleaving peptide) from P, one of the longest yeast promoters with intricate native regulation, in SK1 (also known as KV2695; MATa/α HO gal2 cupS can1R BIO), a diploid yeast strain coexpressing YFP with the FLO11 transcript (56, 57). Gene expression was again bimodal at 38 °C and unimodal at 30 °C (Fig. 2). Overall, these findings suggest that cell fate bifurcation into A cells and R cells at nonoptimal temperatures is a generic, not promoter- or yeast strain-dependent phenomenon.

Considering that the A cell fate also implies protein loss, causing a low-expressing peak and a bimodal distribution, we sought an exception to this effect. In particular, the levels of any R cell fate predictor/marker molecule should be already high even before heating in future R cells. Conversely, the same molecule should be expressed at low levels already before heating in future A cells (Fig. 2). Therefore, the heat-induced loss of an R cell fate predictor protein should not generate a bimodal expression distribution. To test this idea, we investigated a GFP-tagged heat-tolerance marker gene TSL1 (6) at its native locus (Tsl1-GFP S288c; MATa his3Δ1 leu2Δ0 met15Δ0 ura3Δ0) from the yeast-GFP collection. As expected, we found that this yeast strain remained unimodal after heating (Fig. 2). Taken together, these observations support the hypothesis of a cell fate decision between growth arrest and resistance based on molecular marker levels before the onset of heat stress.

Temperature-Dependent Growth Rate of NF Cells.

Knowing that cells undergo a critical decision between growth arrest and resistance, we asked how robust was gene network module function in R cells that continue growing at nonoptimal temperatures. To understand how heating and cooling affect the function of simple network modules in YPH500 R cells, we studied a genomically integrated synthetic biological construct (Fig. 3), the linearizer or NF gene circuit (32). NF consists of the tetracycline repressor (TetR) inhibiting transcription of its own gene and a reporter gene from identical P promoters in an inducer-dependent manner (44, 58). This gene circuit linearizes the dose–response before saturation and reduces the heterogeneity of gene expression compared with similar gene circuits without autoregulation (32). NF can also be a biosensor for molecular effects; for example, deviations of its dose–response from linearity indicates additional feedback (59).

Fig. 3.

Temperature effects on the inducer-dose–response of NF gene expression. (A) Schematic of the negative-feedback (NF) synthetic gene circuit, consisting of the yEGFP::zeoR reporter and the tetR repressor that also regulates its own expression. Repression is relieved by doxycycline (dox). The blunt arrows indicate repression. (B) Experimental (dotted lines) and multiscale Arrhenius model-predicted (solid lines) dose–responses of the high-expressing NF R cell subpopulation mean yEGFP::ZeoR expression (m is a parameter described in , ). Error bars are SEM (n = 3). (C) Dose–response of the coefficient of variation (CV) of R cell yEGFP::ZeoR expression. Error bars are SEM (n = 3). (D) Fluorescence histograms of yEGFP::ZeoR expression in the NF strain at increasing doxycycline concentrations at 12 °C. (E) Fluorescence histograms of the NF strain at 30 °C. (F) Fluorescence histograms of the NF strain at 38 °C. fluor, fluorescence; norm, normalized.

The NF gene circuit is a simple example of a hierarchical network module with regulatory layers through which environmental effects can percolate, making interpretation more difficult compared with constitutive single-gene expression. The NF0 single-gene measurements guided us to first study the exponential growth rates of budding yeast cells with uninduced NF gene circuits (“NF cells” hereafter) at constant temperatures ranging from 5 to 40 °C (). Uninduced NF cells grew optimally at 30 °C and their growth rates decreased similarly to NF0 at other temperatures (). Based on these measurements, we selected 12 °C as the “cold” and 38 °C as the “warm” temperatures to study further, compared with the “standard” 30 °C condition. These temperatures slowed but did not halt cell growth, so that we could still dilute and passage the cell cultures into fresh medium. Earlier evidence that most yeast proteins remain folded below 40 °C (60) and above 5 °C (61) further justified these growth temperature choices.

Armed with an understanding of how temperature affects cell growth with single genes and the uninduced network module, next we studied these effects while altering the inducer level to adjust the strength of regulatory connections. The inducer does not impede yeast growth directly, since doxycycline is harmless to eukaryotes, as confirmed by growing in doxycycline “parental” YPH500 yeast cells devoid of synthetic gene circuit components (). Therefore, to study how temperature interacts with induction to alter NF cell growth, we measured exponential growth rates of NF cells at various inducer (doxycycline) concentrations at warm, cold, and standard temperatures (). We observed that both warm and cold NF cells grew much slower than at 30 °C throughout all induction levels (). To separate the effects of temperature and doxycycline-induced gene expression on growth, we normalized the growth rates in each temperature by their corresponding uninduced values (doxycycline concentration of 0 μg/mL; ). Such normalization eliminates temperature-related effects, so any remaining differences relative to standard temperature should indicate an interaction between induction and temperature. This normalization revealed that heating and cooling increased the slight toxicity of moderate NF induction (). On the other hand, warm conditions abolished the slight toxicity of high NF induction.

Temperature-Dependent Dose–Responses of Noisy Gene Expression in NF Cells.

Next, we asked whether quantitative models of NF cells could predict how temperature affects the dose–response curves of gene expression versus inducer concentrations based on the knowledge we gained so far. To make such predictions, we developed a growth–arrest Arrhenius-type model (, ), which worked well in capturing the temperature-dependent behaviors of the high-expressing NF0 R cells (Fig. 1). We incorporated the growth rates of NF R cells () as well as their arrest rates into a multiscale model, which accounted for population dynamics (arrest and growth of R cells) as well as temperature-dependent intracellular reaction kinetics based on a previously established ordinary differential equation model (32) of the NF gene circuits (Eq. ) (). We first estimated a few unknown parameters by fitting this model to gene expression data at 30 °C and then predicted the dose–response for R cells and thereby for full NF populations at nonoptimal temperatures (, ).

The computational model predicted that cooling should increase the inducer sensitivity of NF dose–responses in R cells relative to the standard temperature and that heating should not affect the gene expression dose–response (Fig. 3). We also studied computationally the temperature dependence of the CV by performing the corresponding stochastic simulations to characterize the widths of the distributions besides their means at various inducer levels (Fig. 4). The models predicted that heating and cooling should leave the NF CV dose–response of high-expressing R cells relatively unaffected (Fig. 4).

Fig. 4.

Revised multiscale models incorporating MD simulation results capture temperature-effects on NF gene circuit function. (A) Probability of TetR DNA-binding domain (DBD) separation distance as a function of temperature from molecular-dynamics (MD) simulations. Distributions of apoTetR with doxycycline bound (−DNA/+dox, denoted −/+) and TetR bound to DNA without doxycycline (+DNA/−dox, denoted +/−). (Inset) Cocrystal structure of TetR bound to DNA (PDB ID 1QPI) with the doxycycline binding pockets shown as pink and gray spheres. The gray ribbons represent TetR (monomers A and B highlighted in red and blue; red and blue spheres are launch points for disordered loops); the sticks colored by element (carbon is cyan) represent DNA. (B) Stochastic dose–response simulations of the yEGFP::ZeoR gene expression mean. m is the temperature-dependent DNA-binding parameter based on MD simulations as described in , . (C) Dose–response of the coefficient of variation (CV) obtained from stochastic simulations of yEGFP::ZeoR expression. Error bars are SEM (n = 3). (D) Stochastically simulated NF gene expression dose–response histograms corresponding to 12 °C. (E) Stochastically simulated NF histograms at 30 °C. (F) Stochastically simulated NF histograms at 38 °C. See for parameters. fluor, fluorescence; norm, normalized; sim, simulations.

Flow cytometry experiments indicated that cooling and heating preserved the NF gene circuit’s mean dose–response linearity (32) nearly up to saturation at all temperatures in R cells (Fig. 3) and at the full-population level (), as predicted computationally. However, cooling increased the slope of linear dose–response and shrank the range of linearity less than heating, in contrast to computational predictions (Figs. 3 and 4 and ). Both heating and cooling lowered fully induced, saturated NF expression (). Neither heating nor cooling altered substantially the CV dose–response of NF R cells (Fig. 3). However, heating elevated the full NF population’s CV dose–response to unusually high levels (), requiring multiscale computational models incorporating the A cells ().

This contrasted sharply all previous studies of the yeast NF gene circuit (32, 59), where gene expression histograms were always unimodal and narrow, maintaining a low CV at all inducer concentrations (Fig. 3). Upon closer inspection, we also found a small low-expressing subpopulation for cold NF cells at high induction (Fig. 3). As for NF0 cells, the CV increase at high temperature stems from gene expression bimodality due to the formation of NF A cells at nonzero induction (Fig. 3). To better understand the process of A cell formation, we used time-lapse microscopy and microfluidics (Movie S1), which again revealed that some heated NF cells randomly arrest growth and concomitantly lose their reporter expression (Fig. 5) along single-cell lineages, similar to NF0 cells. To test A cell viability, we switched the temperature in microfluidic chips to 30 °C but could not directly observe any preheated A cells to resume growth. Nonetheless, we observed small cells percolating through the chip and recovering by initial pseudohyphal growth, within 2 d of restoring the temperature to 30 °C. Accordingly, sorting the lowest-expressing A cells and then batch-culturing them in normal versus heated conditions indicated that some rare A cells do eventually recover at 30 °C, but not at 38 °C (). Overall, heated or cooled cells reach a population-dynamic equilibrium where R cell growth compensates for other R cells converting into nongrowing, nonexpressing A cells, similar to earlier findings (33).

Fig. 5.

Gene expression shutdown and growth arrest in induced NF yeast cells at high temperature. (A) Time courses show the shutdown of GFP reporter expression in three representative induced NF cells at 38 °C and doxycycline concentration of 2 µg/mL. (Inset) Microscopy images of a single NF cell at five different time points corresponding to the time course data indicated by arrows in the main figure. (B) Representative image (bright field and fluorescence overlaid) of induced NF cells trapped in the microfluidic chip for imaging by microscopy at 38 °C and a doxycycline concentration of 2 µg/mL. fluor, fluorescence; norm, normalized.

Additionally, we also tested whether other NF gene circuit variants, incorporating modified regulators and promoters (32), would exhibit analogous temperature-dependent behaviors. Exactly as for NF, cells with these gene circuits had bimodal expression of both the regulator and target reporter genes at 38 °C and unimodal gene expression at 30 °C (), indicating that the cell fate bifurcation did not depend on gene circuit components.

Overall, the mathematical model predicted dose–response changes for cold, but not warm NF gene circuits (Figs. 3 and 4 ). Importantly, the multiscale modeling approach captured the low-expressing subpopulation observed for cooled and heated NF cells (Fig. 4 ) and unimodal expression at optimal temperature (Fig. 4). The remaining disagreement between computational predictions and experimental results suggests some reaction rate-independent (non-Arrhenius) effect of heating. Therefore, we looked for additional effects of heating that could explain this last disparity.

Altered DNA Binding from MD Explains the Functional Effect of Heating NF R Cells.

To elucidate the mechanisms underlying the remaining discrepancy between computational and experimental dose–responses of heated NF R cells, we hypothesized that heating may alter protein structure, affecting transcriptional regulator function. To test this hypothesis, we performed atomistic MD simulations in explicit solvent for the TetR repressor, for which structures are available in the Protein Data Bank ( for details).

TetR functions as a repressor by tightly binding DNA in the major groove, preventing the access of transcriptional machinery. This is possible because the separation of TetR DNA-binding domains (DBDs) is ∼33 Å, which equals the helical pitch of DNA. Previous work (62) has demonstrated that TetR-inducer binding forces the DBDs to move apart, thereby diminishing their propensity to adopt conformations that fit into the major groove of DNA. Therefore, we asked whether temperature might affect the propensity of DNA binding-compatible (or incompatible) TetR DBD conformations.

To address this question, we performed MD simulations on two relevant molecular systems, a TetR:DNA complex (+DNA/−dox; TetR bound to DNA and no doxycycline) and apoTetR (−DNA/+dox; no DNA, but doxycycline is bound to TetR). We simulated these two systems in cold, standard, and warm conditions, testing whether the distributions of apoTetR DBD separation distances (63) change with temperature (). The results (Fig. 4) indicate that without inducer the TetR:DNA complex maintains its dynamics at all three temperatures, sampling from a tight distribution (with small SD) centered at a mean value of 33 Å. This matches the helical pitch of DNA, suggesting that temperature does not alter the DNA-binding affinity of uninduced TetR (63). Likewise, cooled apoTetR adopted a large DBD separation distance, indicating limited DNA binding propensity compared with standard temperature. However, surprisingly, the DBD separation distance distribution of heated apoTetR shifted to lower values. We verified that this shift is not due to DBD unfolding () or doxycycline unbinding (). Unlike for cool and standard temperatures, the warm apoTetR DBD-separation distribution contains the normal TetR:DNA distributions, suggesting that doxycycline-bound apoTetR can abnormally bind to DNA at 38 °C. In summary, the MD simulations suggested that heated apoTetR has an increased probability of binding the operator and repressing gene expression despite the presence of bound inducer.

Motivated by these MD simulation results, we further modified the multiscale mathematical model to also reflect the MD effects, allowing DNA binding of apoTetR in heated R cells. To account for DNA-binding compatible and incompatible conformations of apoTetR, we incorporated temperature-dependent terms into the NF growth–arrest model. This revised model with four temperature effects captured the increased inducer sensitivity for heated NF cells (Figs. 3 and 4). This finding raises new questions about the structural effects of temperature underlying changes in protein–DNA binding and their relevance for gene-regulatory network function.

Multiscale Effects of Heating and Cooling on a PF Gene Network.

To test whether the bimodal gene expression observed at low and high temperatures in induced NF cells (Fig. 3 ) may generalize to other gene networks, we investigated temperature-dependent effects for a second genetic construct, the PF gene circuit (33) ( and Fig. 6). Previously, we employed the PF gene circuit to study how the costs and benefits of gene expression (33) shape population dynamics and thereby genetic evolution (64). The PF gene circuit consists of an rtTA transactivator that activates its own gene and a reporter gene when bound by a tetracycline analog inducer (33).

Fig. 6.

Temperature effects on the inducer-dose–response of PF gene expression. (A) Schematic of the positive-feedback (PF) synthetic gene circuit, consisting of the transactivator rtTA that activates itself and the yEGFP::zeoR reporter when bound by the inducer doxycycline (dox). The pointed arrows indicate activation. (B) Dose–response of the population average yEGFP::ZeoR expression. (C) Dose–response of the overall coefficient of variation (CV) of yEGFP::ZeoR expression. Data points from a subsequent dose–response experiment, performed to characterize the CV peak at 12 °C are shown in cyan. Error bars in B and C are SEM (n = 3). (D) Fluorescence histograms of yEGFP::ZeoR expression for the PF strain at increasing doxycycline concentrations at 12 °C. Inset shows the expression histograms from a subsequent dose–response experiment, which identified equal peaks at doxycycline concentration of 0.06 μg/mL (axes have the same units as the main figure). (E) Fluorescence histograms of PF gene expression at 30 °C. (F) Fluorescence histograms of PF gene expression at 38 °C. fluor, fluorescence; norm, normalized.

As with the NF cells, high and low temperatures reduced PF full-population growth rates in addition to PF induction toxicity (). For PF cells (33), induction-dependent growth rate decline was more pronounced due to activator (rtTA) “squelching” (65) by general transcription factor (TF) sequestration from vital cellular processes (66). Normalizing the growth rates in each temperature as for NF () revealed that low induction still decelerated the growth of warm and cold PF cells, but less than at 30 °C (). This suggests that both heating and cooling might lessen low induction toxicity. On the other hand, at high induction heating seems to worsen induction toxicity, possibly as spurious rtTA interaction affinities toward general TFs increase in an Arrhenius manner.

Heating and cooling appeared to increase the inducer sensitivity of the mean PF dose–response (Fig. 6), which saturated at lower expression levels at both nonoptimal temperatures (), as for the NF cells (Fig. 3 and ). On the other hand, gene expression bimodality emerged at lower induction in cold PF cells, in contrast to warm PF cells, which required higher induction to become bimodal. Considering the cell fate bifurcation of NF and other yeast strains into A cells and R cells, we tested whether these trends were due to the presence of PF A cells in the low-expressing subpopulation. Indeed, there was a low-expressing subpopulation even after exposing heated PF cells to extremely high level of induction with 50 µg/mL doxycycline (), in contrast to the practically unimodal standard PF gene expression already at doxycycline concentration of 4 µg/mL (Fig. 6).

Thus, we calculated the growth rate of R cells from the full PF cell population growth rate () and estimated the fraction of A cells in the low-expressing PF subpopulation at nonoptimal temperatures (). After subtracting the estimated A cell fraction, the heated R cell distributions became more similar, but still not identical to the standard PF distributions, as achieving the same expression pattern still required higher inducer concentrations (). Accordingly, all PF full-population CV dose–responses had a peak whose position and height increased with the temperature (Fig. 6), corresponding to the emergence of bimodal PF histograms (Fig. 6 ). Overall, even apart from the population-dynamic effects of growth arrest, temperature changed the inducer range over which bimodal expression occurred in PF R cells, which is possibly due to Arrhenius effects as well as altered DNA binding dynamics similar to NF R cells.

Discussion

As a first step to mechanistically understand how temperature affects stochastic gene expression in gene-regulatory network modules, we studied the temperature-dependent expression of single genes and synthetic gene circuits (PF and NF). We have identified four different effects of nonoptimal temperatures at various scales: (i) cell fate choice between growth arrest and stress resistance, (ii) diminished growth of the resistant cells, (iii) Arrhenius reaction rates, and (iv) altered protein conformational dynamics. We propose these effects as four pillars of understanding for temperature effects. The second and third of the four pillars are likely most generally applicable, for two main reasons. First, they are based on fundamental and ubiquitous effects: The Arrhenius law applies for all chemical reactions, and growth rate changes affect protein half-lives ubiquitously. Second, they have already been shown to explain temperature-dependent gene network function in bacteria (30) and cell-free systems (31). Nonetheless, it will be interesting to test whether the first pillar (cell fate choice between growth arrest and resistance) might generalize to other stresses and other cell types, considering that findings in yeast often generalize to other organisms (67). The fourth pillar is perhaps the most specific, although it might apply to gene networks incorporating other members of the TetR repressor family (68). Despite sporadic applications of Arrhenius models to investigate how temperature affects genetic oscillators in live bacteria (30) and gene regulation in cell-free systems (69), the generalizability of these models to other networks and other organisms, especially eukaryotes, remained unknown. Moreover, the effects of heating and cooling on stochastic gene expression (noise) have barely been investigated even in prokaryotes (70). Our work on how heating and cooling affects the function of gene circuits controlling stochastic gene expression in yeast fills this double knowledge gap. Nonetheless, future studies will be needed to test the generality of the four pillars across genes, gene networks, and organisms.

Despite the noise-suppressing feature of NF (32, 39, 71), high noise and bimodal expression emerged in the NF gene circuit. Likewise, temperature changes altered the noise of the PF gene circuit, apparently shifting the bimodal region progressively upward in warmer conditions. Explaining these effects required a multiscale modeling framework that combines population dynamics, intracellular biochemical kinetics, and MD. This reveals the multiscale nature of biology and the need for matching models to bridge the experimental and theoretical domains. It also highlights the utility of synthetic biology for developing predictive, quantitative understanding of multiscale biological phenomena. While quantitative understanding of native networks even at a single biological scale can be a struggle, synthetic gene circuits provide new research tools that enable multiscale understanding and prediction of phenomena such as evolution (64) and pattern formation (72). In particular, insights from altered synthetic gene circuit function can reveal structural effects of temperature for low-abundance transcription factor proteins, such as TetR, which challenge even the latest technologies (73). Despite its temperature-dependent slope changes, the NF dose–response remained remarkably linear in R cells, illustrating the robustness that feedback confers to gene network function (74).

Incorporating temperature-dependent kinetic parameters, such as protein synthesis and dilution rates (based on growth rate measurements), enabled us to predict temperature-dependent gene circuit dynamics. Considering the pervasiveness of such temperature effects, these modeling steps should be generally applicable, to predict the temperature-dependent function of other synthetic and even natural gene network modules in various cell types. Considering that the temperature dependence of growth rates is easy to determine, and that protein synthesis rates are directly measurable (75), it will be interesting to explore these possibilities in future, more comprehensive studies.

Determining how TetR-based systems respond to temperature is crucial, for three main reasons. First, they are widely applied to control gene expression in eukaryotic organisms (76). They have been crucially important parts in most synthetic gene circuits developed to date, from the seminal toggle switch (18) and repressilator (19) to logic gates, switches, genetic counters, and other tools. With synthetic biological applications on the horizon, including tens of TetR family members as candidates for novel synthetic biological tools (77), it will be increasingly important to understand how they function at nonoptimal temperatures—as we characterized in the NF and PF strains. Second, TetR is the naming member of a major family of ∼2,353 transcription factor proteins widely present across bacteria, which bear structural and functional similarities (68). Therefore, our findings are potentially relevant for understanding the role, environmental resilience, and functional impact of such proteins within oceanic and soil ecosystems, as well as plant and animal (including human) microbiota exposed to fluctuating temperatures. Third, TetR itself is the key regulator within the transposable TetRA module, the main mediator of tetracycline resistance in bacteria (78), which is temperature sensitive (70). As in the NF gene circuit, the TetR protein transcriptionally represses its own gene and the tetA protein pump, which mediates the most common form of bacterial resistance to tetracycline antibiotics (79, 80). Thus, our findings are generally relevant for the emerging antibiotic resistance crisis, to understand how infectious bacteria carrying TetRA or similar systems resist treatment as temperatures vary (81) due to fever or environmental dissemination.

As synthetic biology moves toward real-world applications, some fundamental challenges will need to be resolved for functional reliability in a diverse range of applications and environments. For example, forward engineering temperature compensation made a synthetic genetic clock robust to temperature changes (30). Inspiration for such designs may come from nature, since temperature compensation appears to be an intrinsic property of natural circadian cycles (82) and certain gene-regulatory network motifs (31). Although temperature-induced changes in synthetic gene circuit function present some challenges, they may also present opportunities to exploit this environmental factor for control purposes. For instance, the period of synthetic genetic oscillators can be tuned by altering temperature (20) in addition to inducer levels and growth media. Our findings suggest that achieving a given gene expression level requires less inducer at low temperatures, which may help to reduce the cost of reagents in scientific and industrial applications. At a fundamental level, the insights gained from these models and experiments open avenues for designing and controlling gene network function and move us toward more fully understanding the dynamics of gene networks in general.

Materials and Methods

Mathematical Modeling.

We developed multiscale models for understanding temperature effects. The models have components at the levels of the cell population, cell interior, and single proteins.

The cell population-level component estimates the growth rate of resistant cells, , and the arrest rate, , based on the dynamics of A and R cell types according to the following ():

The intracellular component investigates the effect of nonoptimal temperatures on gene expression. For NF0 cells, we considered constitutive expression of the yEGFP::ZeoR fluorescent reporter (z) according to the model:

where k(T) and g(T) are temperature-dependent protein synthesis and growth/dilution rates, respectively.

Likewise, to quantitatively understand the temperature dependence of the NF gene circuits dose–responses, we employed previously developed rate equation models. The equations for the NF gene circuit (32) are as follows:

where the variables x, y, r, and z denote the intracellular concentrations of inducer-free repressor (TetR) protein, inducer (doxycycline), inducer-bound TetR protein, and fluorescent reporter (yEGFP::ZeoR) protein. F(x) = θ/(θ + x), where θ is the repression threshold corresponding to an effective repressor-DNA dissociation constant and n is the Hill coefficient. C is a control parameter that describes the rate of doxycycline entry into the cell and is proportional to extracellular inducer concentration. The repressor and reporter-resistance proteins are synthesized at the same rate aF(x) (since both genes are expressed from the same P promoter). Dilution due to temperature-dependent cellular growth of all three variables is g(T); the inducer degrades at a rate f. We neglected the active degradation of TetR and yEGFP::ZeoR due to their relatively long timescales compared with the degradation of the inducer and cell division time in yeast (32). The leaky protein synthesis rate is l and the inducer–repressor binding rate is b.

We added temperature-dependent growth arrest to the NF0 and NF intracellular models by incorporating the experimentally determined low-expressor fraction of cells at each temperature (Fig. 1) into a multiscale version of the Arrhenius model (). This allowed us to investigate the hypothesis that a fraction of the cell population stochastically arrests and shuts off gene expression at high and low temperatures.

The multiscale model with growth arrest accounted for all experimentally observed temperature-dependent results except one: the increased expression sensitivity at high temperature. We rectified this by incorporating temperature-dependent DNA tetR–dox binding (, , and Figs. 3 and 4) that we inferred from submolecular MD simulations.

Based on the results of the MD simulations (, ), we modified the multiscale NF gene circuit model (Eq. ) to incorporate a temperature-dependent parameter m(T) that describes the fraction of apoTetR with a conformational state conducive P promoter binding. This modified model is identical to the model presented in Eq. except F(x) becomes F(x + mr) = θ/[θ + (x + mr)]. Details for a stochastic implementation of this modified model are in .

We obtained the parameters for the 30 °C condition by fitting the NF models to the experimental data using custom Matlab scripts to minimize the sum of squared errors, starting from parameters obtained from (32). We set the dilution rate g for each temperature equal to the mean cellular R cell growth rate obtained from exponential fits to the experimental growth rate data shown in . Note that the Hill function constants have been considered temperature independent, as they are ratios of rate constants and are specific for a certain reaction (83). The Arrhenius scaling of reaction rates for the 12 and 38 °C temperature conditions is discussed in .

For details on MD simulations, see .

Strains and Media.

The haploid Saccharomyces cerevisiae strain “YPH500” (MATα ura3-52 lys2-801_amber ade2-101_ochre trp1-Δ63 his3-Δ200 leu2-Δ1; congenic with the standard laboratory strain S288c) was used as a model organism throughout most of this study. “NF0” (YDN-G1GZmbh), “NF” (YDN-G1GZmbh-G1Tbt), and “PF” (YDN-T2dGZmxh-T2dMFot) cultures were grown in synthetic drop-out (SD) medium with the appropriate supplements (NF0, SD −his +ade; NF and PF, SD −his −trp +ade) to maintain selection and supplemented with sugars (2% glucose or 2% galactose) as described below (all reagents from Sigma). Doxycycline (Fisher Scientific) was used to induce the PF and NF gene circuits. Additionally, experiments were performed on “2-Color NF” (YDN-G1GZmbh-G1TCbt), “2-Color NF-T123” (YDN-T1GZmbh-T1TCb), “TBR1” (YLC-TBR1-G1Gh), “Kv2695” (Sk1 FLO11-2A-YFP; K. J. Verstrepen laboratory, KU Leuven, Leuven, Belgium), and “TSL1” (Tsl1-GFP S288c; S. F. Levy laboratory, Stanford University, Palo Alto, CA) yeast strains; the former three strains were grown in SD media (2-Color NF, SD −his −trp +ade; 2-Color NF-T123 and TBR1, SD −his +ade) and the latter two strains grown in YPD and YPGal media. YPH500 (Stratagene), the genetic ancestor to the NF0, NF, PF, 2-Color NF, 2-Color NF-T123 strains, was also cultured in YPD and YPGal media.

Cell Culture.

Well-isolated single yeast colonies were picked from plates and incubated overnight in synthetic drop-out medium supplemented with 2% glucose at 30 °C. Twelve hours later, the 1-mL cell suspensions were diluted to the concentration 1 × 106 cells per mL—concentration estimated using a Nexcelom Cellometer Vision cell counter (Nexcelom Bioscience)—in fresh medium supplemented with 2% galactose. Triplicate cultures were grown for 48 h in galactose medium with the appropriate concentration of doxycycline (Acros Organics) for each condition to allow gene expression levels to stabilize (32, 33). Cell density was measured (Nexcelom Cellometer Vision) and resuspended to a concentration 1 × 106 cells per mL every 12 or 24 h depending on the experimental condition to keep them in log growth phase.

Flow Cytometry and Fluorescence Microscopy.

Population-level gene expression was read on the BD Accuri flow cytometer (Becton Dickinson) after 48 h of doxycycline induction. We considered gene expression being stable when GFP distributions for 30 and 38 °C experiments did not change between 24- and 48-h measurements (NF cells: Fig. 3 and ; PF cells: Fig. 6 and ), and similarly, the distributions for 12 °C did not change between 48- and 72-h measurements (NF cells: Fig. 3 and ; PF cells: Fig. 4 and ).

For fluorescence microscopy, cells were cultured for 2 d at 38 °C before imaging () and then transferred to a modified high-throughput yeast aging analysis microfluidics chip (49). Images were acquired on a Nikon computer-controlled motorized inverted fluorescence TIE Eclipse microscope (Nikon) equipped with a Nikon DS-Qi2 camera. Composite images were prepared in NIS Elements (Nikon).

Growth Rate Measurements.

Population growth rate between subsequent resuspensions was estimated by linear fits to log-transformed cell counts (inferred from cell density measurements, described in , and culture volume).

Data Processing and Analysis.

Raw flow cytometry data files were read into Matlab (Mathworks) using the Matlab script fca_readfcs (Matlab Central) for plotting and analysis. A small gate was applied to the forward-scatter and side-scatter data to minimize the contribution of extrinsic noise due to cell cycle phase, cell size, and age (84), and exclude doublets, dead cells, and cellular debris from the analysis. To eliminate small numbers of mutated cells that may have lost the integrated construct (due to homologous recombination) or rare cells left over from previous samples (not eliminated by flow cytometer), cells with log fluorescence deviating more than 3 SDs from the mean were considered outliers and discarded from the analysis (32, 33, 59). Time-lapse microfluidics images were analyzed in Matlab using custom scripts. All data and Matlab scripts are available at https://openwetware.org/wiki/CHIP:Data.