Genetic Mutations That Drive Evolutionary Rescue to Lethal Temperature in Escherichia coli.

Evolutionary rescue occurs when adaptation restores population growth against a lethal stressor. Here, we studied evolutionary rescue by conducting experiments with Escherichia coli at the lethal temperature of 43.0 °C, to determine the adaptive mutations that drive rescue and to investigate their effects on fitness and gene expression. From hundreds of populations, we observed that ∼9% were rescued by genetic adaptations. We sequenced 26 populations and identified 29 distinct mutations. Of these populations, 21 had a mutation in the hslVU or rpoBC operon, suggesting that mutations in either operon could drive rescue. We isolated seven strains of E. coli carrying a putative rescue mutation in either the hslVU or rpoBC operon to investigate the mutations' effects. The single rescue mutations increased E. coli's relative fitness by an average of 24% at 42.2 °C, but they decreased fitness by 3% at 37.0 °C, illustrating that antagonistic pleiotropy likely affected the establishment of rescue in our system. Gene expression analysis revealed only 40 genes were upregulated across all seven mutations, and these were enriched for functions in translational and flagellar production. As with previous experiments with high temperature adaptation, the rescue mutations tended to restore gene expression toward the unstressed state, but they also caused a higher proportion of novel gene expression patterns. Overall, we find that rescue is infrequent, that it is facilitated by a limited number of mutational targets, and that rescue mutations may have qualitatively different effects than mutations that arise from evolution to nonlethal stressors.

Significance

In lethal environmental conditions, population survival is thought to occur only through ecological means (e.g. migration). However, populations can respond genetically and evade extinction in a process known as evolutionary rescue. Although rescue has been documented previously, the genetic mutations that drive rescue and their consequences have rarely been thoroughly studied. Here, we performed experiments to identify the adaptive mutations that allow Escherichia coli populations to survive high temperature beyond its thermal niche.We identified two adaptive pathways that facilitate population rescue,measured their fitness effects and assessed their effects on gene expression. Ourwork suggests that rescuemutations may have different effects than mutations arising from evolution in nonlethal conditions.

Introduction

Under severe environmental stress, a population will decline rapidly and may face extinction. However, populations can adapt genetically; if an individual appears with an adaptation to the severe stress, the population may recover. This process of decline and recovery results in a U-shaped pattern of population dynamics that defines the phenomenon of evolutionary rescue (Bell 2017). It is important to understand the frequency and dynamics of rescue events, both because they affect our understanding of species’ survival and also because they have practical implications for medicine, agriculture, and conservation biology. For example, evolutionary rescue drives some of the dynamics of bacterial antibiotic resistance. When exposed to potentially lethal concentrations of antibiotics, a bacterial population declines, but a resistance mutation can restore population growth (Orr and Unckless 2008; Baquero and Cantón 2017). This evolutionary response often occurs because antibiotics target a specific enzyme or structure (e.g. the ribosome), so that a single beneficial mutation inhibits the antibiotic’s mechanism of action (Blair et al. 2015). Similar dynamics contribute to fungicide and pesticide resistance in agriculture, where there is often a simple genetic basis to resistance (Délye et al. 2013; Lucas et al. 2015).

Sometimes environmental challenges affect more complex physiological traits, and this issue has been addressed to some extent in experimental studies of evolutionary rescue. In one study, yeast populations were grown in a lethal concentration of salt (Bell and Gonzalez 2009). Salt tolerance in yeast is a complex, polygenic trait (Dhar et al. 2011), suggesting that evolutionary rescue could result from mutations in one of several genes or perhaps even require multiple genic changes. Studies have shown that adaptive mutations can rescue yeast populations from lethal salt conditions, but such rescue occurs infrequently, and the probability of rescue varies by population size (Bell and Gonzalez 2009). Similarly, green alga (Chlamydomonas) has also been used to study evolutionary rescue to another complex trait, low-light conditions (Bell 2013), showing again that rescue is infrequent and depends on population characteristics. Together, these studies highlight that rescue can alter evolutionary outcomes for complex traits. However, neither set of studies identified the genetic basis of rescue (Bell and Gonzalez 2009; Bell 2013; Gonzalez and Bell 2013), which is an important precursor for understanding the dynamics of rescue and its underlying mechanisms.

Here, we study evolutionary rescue in Escherichia coli that has been challenged with a lethal temperature. Temperature is a complex environmental variable because it governs the rates of biological reactions that underlie respiration, growth, and reproduction (Somero 1978; Cooper et al. 2001). Furthermore, characterizing the evolutionary response to severe thermal stress is important for understanding adaptation to global climate change (Holt 1990). Previous work has investigated adaptation to both nonlethal and lethal heat stress using E. coli. As an example of nonlethal stress, Tenaillon subjected E. coli to a high but sustainable temperature (42.2 °C) and found >1,000 putatively adaptive mutations, illustrating the genetic diversity of adaptive responses. These mutations occurred within dozens of genes, but there were also clear patterns. Mutations were especially frequent in genes that modify transcription, such as the RNA polymerase (RNAP) subunit β (rpoB)gene. Some of these rpoB mutations conveyed a fitness benefit in high heat but fitness trade-offs at lower temperatures (<20.0 °C), indicating antagonistic pleiotropy (Rodríguez-Verdugo et al. 2014).

Previous work has also subjected E. coli to lethal temperatures that resulted in rare rescue dynamics (Bennett and Lenski 1993). For example, Mongold characterized patterns of evolutionary recovery at 44.0 °C, using E. coli strains that had been adapted to 32.0, 37.0, and 41.0–42.0 °C (Mongold et al. 1999). They found that rescue events at 44.0 °C occurred in 8% of populations but only in populations derived from ancestors that had been previously adapted to high temperature (41.0–42.0 °C), suggesting that preadaptation contributes to evolutionary rescue. Moreover, they found that the rescued populations exhibited a fitness cost at elevated, but nonlethal temperatures, suggesting that at least some rescue mutations were antagonistically pleiotropic. Here again, however, the underlying adaptive mutations were not identified.

In this study, we perform E. coli growth experiments to better understand the dynamics, mechanism, and fitness consequences of evolutionary rescue. Beginning with an ancestor derived from a single colony, we carry out replicated evolution experiments at 43.0 °C, which typically results in population extinction under our growth conditions. After observing and noting the frequency of rescue events, we identify mutations within the rescue populations. With these mutations in hand, we ask the following three sets of questions. First, can these mutations drive evolutionary rescue to lethal temperature in E. coli? Discriminating between driving and hitchhiking mutations is a major challenge in evolutionary biology (Rosenzweig and Sherlock 2014), and hence unambiguous identification of drivers is an important goal. Second, what is the fitness effect of driver mutations, and do they have trade-offs that affect their population dynamics? We are specifically interested in antagonistic pleiotropy, a phenomenon that has been shown to be common (Williams 1957; Cooper and Lenski 2000; MacLean et al. 2004) but not universal in evolution experiments. Finally, can we glean any insights into the molecular effects and mechanisms of rescue? To do so, we study gene expression changes introduced by driver mutations, to try to better understand their downstream effects. We also assess whether gene expression shifts back toward an unstressed physiological state (Carroll and Marx 2013) or toward novel expression patterns.

Materials and Methods

Evolutionary Rescue Experiments

A frozen glycerol stock was prepared from a single colony of E. coli B strain REL1206 possessing a neutral Ara− marker. This strain had been propagated previously at 37.0 °C for 2,000 generations in Davis minimal medium supplemented with glucose at 25 mg/l (DM25) and was thus adapted to the growth medium (Lenski et al. 1991). To isolate the single colony, REL1206 was streaked from frozen onto a tetrazolium–arabinose (TA) plate and incubated overnight at 37.0 °C. The single colony was inoculated into Luria-Bertani medium (LB) and grown overnight. To prepare a frozen reference stock, 900 μl of culture was mixed with 900 μl of 80% glycerol and frozen at −80 °C. We term this REL1206 frozen stock the “rescue ancestor” (fig. 1). A backup rescue ancestor stock was prepared from the same LB culture.

. 1.

Evolutionary rescue experimental design and dynamics. (A) Experimental design for producing and observing rescue events. Bacteria were propagated from frozen through two flasks to acclimate them and to produce enough cells for experimental replication. Samples of flask culture were transferred to 44 replicates that were propagated through 1:100 serial dilution for 5 days. This procedure was repeated across 7 different weeks. (B) Population cell densities over time. Most populations went extinct over the course of 5 days. A total of 26 rescue events were observed across the 3rd, 4th, and 5th days of growth. The timing of rescue events was determined by the day at which cell density increased by an order of magnitude over the previous day. Populations possessing rpoBC mutations are indicated by rectangles. The three populations possessing duplications are circled. (C) Genome-wide distribution of mutations in rescue populations. Populations 1–26 are labeled on the left and different weeks are separated into groups and labeled at the right. Mutations are colored by their frequency in the population according to the scale at the right. Synonymous, nonsynonymous, indel, and intergenic mutations are represented by squares, circles, triangles, and diamonds, respectively. Only mutations at frequencies >10% are shown. Mutations occurring in more than two populations are labeled at the top.

As is common practice (Bennett and Lenski 1993; Lenski and Travisano 1994; Rodríguez-Verdugo et al. 2014; Hug and Gaut 2015), we first acclimated the rescue ancestor (REL1206) to mild laboratory conditions to allow it to recover from being frozen. The rescue ancestral stock was inoculated into 100 ml LB and grown for 8 h in an Infors HT Minitron incubator at 37.0 °C and 120 rpm (fig. 1). 10 μl of this culture was then inoculated into 100 ml DM25 and grown for 24 h in an Infors HT Minitron incubator at 37.0 °C and 120 rpm. We inoculated 100 μl of the 37.0 °C, DM25 culture into each of 44 culture tubes containing 9.9 ml DM25. An additional four culture tubes containing 9.9 ml DM25 were used as contamination controls and cell density blanks. The total set of 44 tubes was placed into an Innova 3100 water bath shaker (New Brunswick Scientific) and grown for 24 h at 120 rpm and at the experimental temperature of 43.0 °C. Tube cultures were serially propagated over the course of 5 days by inoculating 100 μl of culture into 9.9 ml DM25 after each day of growth (fig. 1). Note that the 44 inoculated tubes did not constitute independent experiments, because they all derived from the same overnight, 100 ml DM25 culture. However, the procedure was repeated independently across 7 weeks, for a total of (7 weeks × 44 =) 308 populations. Our interpretations make use of both the independent and nonindependent features of the design.

To measure cell densities from individual populations, 50 μl of each culture (including blanks) were inoculated into cuvettes containing 9.9 ml Isoton II Diluent (Beckman Coulter) on each of the 5 days. These samples were analyzed using a Multisizer 3 Coulter counter (Beckman Coulter) to determine cell densities (particles/ml). The average particle density of the four blanks was subtracted from each sample’s cell density. Of 308 experimental populations, 26 were excluded due to technical failure with the Coulter counter, leaving 296 measured populations. We defined a rescue event as a population whose cell density increased by at least an order of magnitude over the previous day’s measurements. Following cell density measurement, rescue populations were saved as frozen glycerol stocks on their initial day of recovery, as well as any subsequent days. To prepare frozen stocks of rescue recovery events, 900 μl of culture was mixed with 900 μl of 80% glycerol and frozen at −80 °C.

DNA Extraction and Population Sequencing

Most samples for DNA extraction and sequencing were derived from day 5 of the experiment, but two samples (populations No. 1 and No. 19) were derived from day 4. Each of the frozen rescue populations was inoculated into ten separate culture tubes containing 9.9 ml DM25 and incubated in an Innova 3100 water bath shaker (New Brunswick Scientific) overnight at 37.0 °C and 120 rpm. Cells from all ten tubes were pooled, and genomic DNA was extracted from these samples using Wizard Genomic DNA Purification Kits (Promega). Pooling was employed to filter new mutations that might have risen during the process of recovery. DNA from the rescue ancestor was extracted in the same manner but pooled from four tubes rather than 10. Genomic DNA libraries were prepared using the TruSeq DNA PCR-Free Library Preparation Kit (Illumina). The 26 rescue populations were multiplexed and sequenced in two lanes of an Illumina HiSeq 2500 in rapid mode at UC Irvine’s Genomics High-Throughput Facility. Two ancestral samples—one working stock (fig. 1) and one backup stock—were also sequenced on an Illumina HiSeq 3000 at the Bioinformatics Core Facility at the UC Davis Genome Center.

Mutations and mutation frequencies were called using breseq (Deatherage and Barrick 2014) in polymorphism mode, using the E. coli B REL606 genome as a reference, which differs from REL1206 in six positions that were excluded from our analyses (Barrick et al. 2009; Tenaillon et al. 2012). In theory, breseq provides information about duplications and deletions by reporting novel junctions. No evidence of novel junctions or sequencing coverage was found for large deletions in our data set, but breseq did provide some novel junction evidence for the presence of large duplications. To assess duplications more formally, we compared unique reads (mapping quality >5 in samtools 1.3) across 10 kb regions of the genome, defining duplications as regions with more than twice the average genome coverage.

Isolating Single Mutants from Rescue Populations

To purify isolates carrying a single mutation in either the hslVU or rpoBC operon, we selected populations that had a single fixed mutation in hslVU or rpoBC; that is, occurring at a frequency of 85% or greater as called through the breseq analysis. These populations were streaked from frozen stock onto TA plates and incubated at 37.0 °C. Multiple single colonies were picked per line and subsequently purified on new TA plates. Purified isolates were then grown in LB and incubated at 37.0 °C and 120 rpm. This culture was used to prepare frozen glycerol stocks of the purified isolates.

For each purified isolate, we designed PCR primers (https://www.idtdna.com; last accessed May 2018) for all mutations found in the population at >10% frequency. We submitted the PCR products for Sanger sequencing to determine mutation presence or absence for each screened gene. From the PCR and Sanger sequencing results, we identified putative single mutant genotypes. Isolates that were positive for any nonfixed (background) mutations were eliminated from further study. To determine the validity of these putative single mutant genotypes, we performed whole-genome sequencing. Total genomic DNA from the putative single mutants was extracted using the Promega Wizard Genomic DNA Purification kit, and DNA concentrations were measured with Qubit dsDNA HS Assay kits. Genomic DNA libraries were prepared for sequencing using the Illumina Nextera DNA Flex Library Preparation kit. Libraries were multiplexed and sequencing was carried out on a single lane of Illumina HiSeq 4000 at the UCI Genomics High-Throughput Facility. To call mutations in these isolates, the Illumina reads were mapped against the reference sequence, REL606, using breseq as described above (Deatherage and Barrick 2014).

Relative Fitness Measures

To measure relative fitness, we competed populations or single mutants against the ancestral line by growing them together in the same culture tube. We used REL1207 as the ancestral strain, because it is identical to the REL1206 strain except for a neutral marker, Ara+. High temperature fitness assays were performed at 42.2 °C, because REL1207 does not survive at higher temperatures in our system.

To perform assays, we revived REL1207 and either the rescue populations or single mutants from frozen into 10 ml LB and incubated at 37.0 °C with 120 rpm. The next day, we diluted the overnight cultures 100-fold in saline and transferred 100 μl of this dilution into 9.9 ml DM25 media. This was then incubated at 37.0 °C with 120 rpm to acclimate from frozen conditions (Bennett and Lenski 1993). After 24 h of incubation, the cultures were transferred to fresh DM25 media and incubated at 42.2 °C to acclimate to high temperature stress. The following day, we mixed the ancestral strain and a rescue mutant or population 9:1 into sterile DM25 media. This mixture was plated onto TA solid media to count the initial cell densities before competition. This mixture was incubated at 42.2 °C with 120 rpm. The cells were left to compete for 24 h, and we quantified the final cell densities of the ancestor and rescue line by plating on TA plates and counting colonies. To perform competitions at 37.0 °C, we began the competition on the day following the first acclimation step in DM25 and mixed the ancestral and rescue lines at a 1:1 ratio.

To quantify relative fitness, wr, we used the methods as in Lenski and Tenaillon . The fitness of a mutant or population relative to the ancestor is estimated as: , where and represent the initial cell densities of the mutant (or population) and the ancestor before competition, and and represent the final cell densities after 1 day of competition.

RNA Harvest, Isolation, and Sequencing

To harvest cells for RNA extraction, we grew the single mutants to the mid-exponential phase of their growth curve. To do so, we acclimated the single mutants from frozen stock in 10 ml LB media at 37.0 °C with 120 rpm in an Innova 3100 water bath. We then diluted these cultures 10,000-fold into DM25 media and incubated the cultures at 37.0 °C with 120 rpm. Following 24 h of incubation, the cultures were diluted 1,000-fold into DM25 media and acclimated to either 42.2 °C or 43.0 °C with 120 rpm for 24 h as is customary to acclimate cells to stressful temperatures (Bennett and Lenski 1993; Rodríguez-Verdugo et al. 2014). The following day, the growth curve was started by transferring 100 μl of the culture to 24 tubes with 9.9 ml of DM25 and incubated at either 42.2 °C or 43.0 °C. Cell density was measured using the Multisizer 3 Coulter counter (Beckman Coulter) in volumetric mode by diluting 50 μl of cell culture into 9.9 ml of Isoton II diluent. We measured the cell density every 30 min following the first 5 h of growth until the cells reached the mid-exponential growth phase based on the electronic counts. Cells were concentrated through vacuum filtration of 150–200 ml of culture onto cellulose nitrate membrane filters with 0.2 μm pore size. The cells were washed off from the filters and pelleted for storage at −20.0 °C in a mixture of 2 ml of Qiagen RNA Protect Bacterial reagent and DM25 media. Three replicates per single mutant line were harvested for both temperatures (7 mutants × 2 temperatures × 3 replications = 42 samples), and three replicates of the ancestral line, REL1206, were harvested at 42.2 °C, for a total of 45 RNAseq samples.

The cell pellets were thawed and treated with lysozyme for 5 min before extracting total RNA using Qiagen RNeasy kits. RNA concentrations were measured with Qubit RNA HS assay kits and RNA quality was assessed by running an Agilent RNA-Nano chip on a bioanalyzer. We enriched for mRNA by the removal of rRNA using NEBNext rRNA depletion kits for bacteria. We prepared the RNA for Illumina sequencing using the NEB Ultra II Directional RNA Library Prep kit. All samples were uniquely barcoded and multiplexed for sequencing with Illumina NovaSeq at the UCI Genomics High Throughput Facility.

Gene Expression Analyses

RNA sequencing reads from our study and previously sequenced reads from REL1206 grown at 42.2 and 37 °C from Rodríguez-Verdugo et al. (2016) were filtered with a custom Perl script to a quality cut-off of 20. The filtered reads were then mapped to the REL606 reference sequence using BWA version 0.7.8 with default parameters (Li and Durbin 2009). Uniquely mapping reads were used as input into HTSeq, which counts the number of uniquely mapped reads to annotation features (Anders et al. 2015). Analysis of the RNAseq counts was carried out in R (R Core Team 2019). We normalized the RNAseq counts and identified differentially expressed genes (DEGs) using the DESeq2 package (Love et al. 2014). We followed previous studies (Rodríguez-Verdugo et al. 2016; González-González et al. 2017) by identifying DEGs as significant at Padj < 0.001 and also exhibiting log2-fold change >2 between samples. Gene Ontology (GO) enrichment analyses were performed at the online website (http://geneontology.org; last accessed July 2020) using E. coli as the reference list (Ashburner et al. 2000; The Gene Ontology Consortium 2019).

Because we used new and previously published RNAseq data for expression analyses, we were concerned about the potential for batch effects. To assess batch effects, we compared analyses comparing single mutants at 42.2 °C to previous REL1206 data at 42.2 °C (n = 2; Rodríguez-Verdugo et al. 2016) and to our new RNAseq data (n = 3) of REL1206 at 42.2 °C. The new data resulted in 20% more detected DEGs than the old data, perhaps reflecting differences in power with different sample sizes (n = 3 vs. 2). Importantly, however, 93% of DEGs were shared between the two analyses, and the two data sets led to qualitatively identical GO-enrichment analyses. Based on this comparison, we concluded that batch effects did not dramatically alter overall conclusions about the types and direction of genic shifts in expression. We therefore combined old and new REL1206 42.2 °C samples, so that all reported comparisons with the 42.2 °C ancestor were based on n = 5 replicates.

Once detected, changes in DEG gene expression were categorized into one of four directions (restored, reinforced, novel, or unrestored) as previously described (Carroll and Marx 2013; Rodríguez-Verdugo et al. 2016; González-González et al. 2017). These directions represent the change in gene expression of a rescue mutant relative to the ancestor’s gene expression at 42.2 and 37.0 °C, where 42.2 °C represents a stressed state for the ancestor and 37.0 °C is an unstressed state (supplementary table S1, Supplementary Material online). Briefly, a gene was restored if the ancestral expression level was significantly different from itself at 37.0 and 42.2 °C, and the mutant expression level was significantly different and in the opposite direction to that of the ancestral gene expression at 42.2 °C. A gene was reinforced if the ancestral expression level whereas stressed at 42.2 °C was significantly different from its unstressed expression level at 37.0 °C, and the mutant’s expression level of the gene was significantly different and exaggerated in the same direction to that of the ancestral gene expression at 42.2 °C. A gene had novel expression if the ancestral expression level was not significantly different from itself at 37.0 and 42.2 °C, but the mutant had significantly differential expression to the ancestral expression level at both temperatures. Finally, a gene was unrestored if the ancestral expression level was significantly different from itself at 37.0 and 42.2 °C, and the mutant did not have a significant difference in expression to the ancestral gene expression level at 42.2 °C.

Results

Mutations Associated with Rescue Events

We ran experiments that started with an overnight culture of the REL1206 ancestor at 37.0 °C in low nutrient DM25 media (fig. 1). We then transferred 100 μl of the overnight culture into 44 tubes of fresh DM25 media. These 44 cultures represented distinct populations which we maintained for 5 days by 1:100 daily serial dilution at 43.0 °C in a precisely controlled shaking water bath (fig. 1). This experiment was repeated over 7 separate weeks, for a total of 308 (=7 × 44) experimental populations. Of these, 296 populations were monitored for cell density over a period of 5 days to determine whether the population went extinct or rebounded to rescue. Altogether, we identified rescue events in 6 of the 7 weeks and 26 populations. Thus, the frequency of rescue was 8.8% of populations (i.e. 26 of 296). Among the rescue events, 3 were detectable on day 3, 12 more on day 4, and the rest on day 5 (fig. 1). Neither the number nor the timing of rescue events were correlated with initial cell densities (r = −0.24, P = 0.61; r = 0.50, P = 0.32), suggesting that results were not driven by variation in initial conditions across weeks.

To characterize the genomic changes associated with rescue events, we sequenced each of the 26 rescued populations and identified the frequencies of mutations. We focused on mutations that reached near-fixation, which we defined as >85% frequency. The 26 populations contained 1.8 fixed mutations on average, but 11 populations had just one fixed mutation, making these genetic changes the likely drivers of population recovery. One rescue population (No. 2) evolved a mutator phenotype due to a small deletion in the mutT gene and contained six fixed mutations, the most of any population in the experiment (fig. 1). Four populations contained two distinct large duplications based on their sequencing coverage profiles; one duplication included the groEL and groES genic region, and the second contained a duplication that included the hslVU operon.

In total, the 26 populations yielded 29 distinct point or small indel mutations within 20 different genic and intergenic regions (table 1). However, three regions were particularly notable. The first was the clpA/serW intergenic region, in which the same point mutation appeared across 4 separate weeks (table 1). Note, however, that this mutation always appeared with other fixed mutations, making it unclear whether it was sufficient to drive rescue. The second was the rpoBC operon, where four distinct mutations were identified across 5 of 7 weeks and 6 of 26 populations. All of these point mutations caused nonsynonymous changes, including one rpoC mutation (W1020G) that was fixed across 3 separate weeks. Two of the rpoBC mutations (rpoB H447L, rpoC W1020G) were the only fixed mutations in at least one population, suggesting that they were sufficient to rescue a population.

Table 1

Fixed Mutations in Rescued Populations

Gene(s)a	Week–Populationb	Positionc	Mutationd	Mutation Typee
arcB	1–2	3285997	T-->G	E755A (GAA-->GCA)
clpA/serW	1–2, 2–6, 2–7, 3–26, 3–16, 5–20	942604	A-->G	Intergenic (+359/+339)
ECB_00530	3–16, 5–21, 6–24	573246	G-->C	L81L (CTC-->CTG)
ECB_00530	2–5	573229	T-->G	N87T (AAC-->ACC)
ECB_02812/ECB_02813	1–2	3012076	T-->G	Intergenic (−516/−58)
hslU	3–9, 3–18, 3–19	4100115	+C	Coding (1116/1332 nt)
hslU	2–5, 2–8	4100512	Δ1 bp	Coding (719/1332 nt)
hslU	3–11	4101052	C-->T	G60D (GGT-->GAT)
hslU	2–6, 2–7	4101159	C-->A	K24N (AAG-->AAT)
hslU	3–14	4100743	A-->C	L163R (CTG-->CGG)
hslV	2–4, 5–21	4101760	+T	Coding (11/531 nt)
hslV	3–17	4101363	Δ2 bp	Coding (408/531 nt and 409/531nt)
hslV	3–12, 3–13	4101568	T-->G	H68P (CAT-->CCT)
hslV/ftsN	1–1, 6–25	4101844	Δ1 bp	Intergenic (−74/+19)
insE-1/serX	2–4, 2–7	1111967	A-->G	Intergenic (−236/+166)
lnt	2–4	671609	A-->C	G456G (GGT-->GGG)
mrdA	1–2	649901	T-->G	I301L (ATC-->CTC)
mutT	1–2	114029	Δ1 bp	Coding (182/390 nt)
pepA	2–8	4468692	T-->G	T163P (ACC-->CCC)
rhsE	3–18, 6–24	1500351	T-->G	G25G (GGT-->GGG)
rpoB	1–3	4162195	A-->T	H447L (CAC-->CTC)
rpoB	1–2	4163133	A-->C	N760H (AAC-->CAC)
rpoC	5–20	4165883	A-->G	D308G (GAT-->GGT)
rpoC	3–16, 6–24, 7–26	4168018	T-->G	W1020G (TGG-->GGG)
rtcA	2–4	3484800	G-->A	S217F (TCC-->TTC)
secF	3–10	398683	+GGT	Coding (756/972 nt)
ybgG/cydA	2–8	751779	T-->G	Intergenic (+286/−561)
yghS	2–8	3067389	T-->G	H219P (CAC-->CCC)
ynfL	2–4	1646703	T-->A	I29F (ATT-->TTT)

Gene names as defined in the REL606 annotation. Bolded names represent single mutations isolated for further study (see table 2).

Provides information about the week and population in which the mutation was fixed. This column shows, for example, that the same clpA/serW intergenic point mutation was present across four independent weeks and six total populations.

Location of mutations in the REL606 reference.

Provides mutation type from REL1206 for point mutations. + and Δ represent an insertion and deletion relative to REL1206.

Provides information about codon change for nonsynonymous mutations. Coding variants with + and Δ represent frameshifts.

Table 2

Seven Single Mutants and Their Relative Fitness Values

	Relative Fitness 42.2 °C		Relative Fitness 37.0 °C
Mutant	w r	w r P value	w r	w r P value
hslU frameshift	1.317	0.003	0.933	0.033
hslU G60D	1.060	0.073	0.993	0.378
hslU L163R	1.419	4.88E−05	0.992	0.272
hslV frameshift	1.115	0.0004	0.967	0.062
hslV H68P	1.192	0.011	0.966	0.195
rpoB H447L	1.303	0.0054	0.905	0.006
rpoC W1020G	1.268	0.0002	0.993	0.257

Finally, the most mutations were observed in the hslVU operon, which encodes a heat shock protease system (Missiakas et al. 1996; Bochtler et al. 2000). In addition to the duplication of this region mentioned previously, for which both copies apparently had a 1 bp indel frameshift (table 1, populations 5–21), the operon had eight distinct point or indel mutations across 5 of the 7 weeks. Another fixed mutation altered the 3′ intergenic region of this operon (table 1). Altogether, hslVU mutations were found in 62% (16/26) of rescued populations and were the only fixed mutation in at least seven populations. Five of the mutations within hslVU caused frameshifts, suggesting that interruption of function was adaptive. Interestingly, the populations with fixed rpoBC mutations were distinct from those with hslVU mutations; no populations contained fixed mutations in both operons, even though mutations in both operons were identified during weeks 3 and 5 (table 1 and fig. 1).

We evaluated four additional features of fixed mutations, focusing on populations that contained potential driver mutations in rpoBC and hslVU. First, because rescue events occurred at different times during the course of the 5-day experiments, we assessed whether the identities of fixed mutations were related to the day of rescue. We found no obvious relationship (Mann–Whitney U test comparing hslVU and rpoBC populations: P = 0.34). Second, we tested for a relationship between fixed mutations and cell densities after population recovery; rpoBC populations had significantly higher cell densities than all other recovered populations on both days 4 and 5 (Mann–Whitney U test: day 4 P = 0.0044, day 5 P = 3.71 × 10−4) (fig. 1). Third, to assess whether the fixed mutations arose as a consequence of thermal stress, we also sequenced two control cultures that were maintained at 37.0 °C in DM25 (see Materials and Methods section). Both controls were sequenced to >2,000×, but we found no fixed mutations relative to the REL1206 genome, only three mutations present at >10% frequency, and an average variant frequency of 1.3%. Of the three mutations at >10%, two were not shared with any of the rescue populations. The remaining mutation was a four-nucleotide indel within the ECB_01992 gene which is likely to be hypermutable because it is part of a motif of seven four-nucleotide repeats (Tenaillon et al. 2016). The four-nucleotide indel was found at frequencies of 19% and 21% in the two ancestral samples and also between 32% and 34% frequency across 8 of the 26 rescue populations (supplementary table S2, Supplementary Material online). The fact that it was not fixed in any population suggests it was likely not adaptive.

Fitness Properties of Single Mutations in hslVU and rpoBC

Repeated mutations in rpoBC and hslVU across weeks and populations suggested that specific mutations in these operons drive evolutionary rescue. To verify this conjecture, we isolated seven clones containing single mutations in either the rpoBC operon or the hslVU operon (tables 1 and 2).

For each of the seven single mutants, we measured their relative fitness (wr) against the ancestor using competition assays. To measure wr we competed the ancestor against each of the single mutants at 42.2 °C, because the ancestor does not survive at higher temperatures. Six of the seven mutants conferred a significant advantage (wr > 1.0; fig. 2 and table 2), and collectively they had a 24% average wr increase, (one-tailed t-test: P = 3.44 × 10−12). The two rpoBC mutants had 26% and 30% fitness advantages, whereas the hslVU mutants ranged from an estimated 5% to 41% advantage. We note that wr increases for single mutants were nearly identical to those based on wr estimates based on competing population samples against the ancestor (supplementary table S2, Supplementary Material online). On average, the 26 rescued populations had a wr increase in 24% relative to the ancestor, with an average of 28% (n = 7; range: 14–35%) and 23% (n = 16; range: 4–47%) for populations that contained fixed rpoBC and hslVU populations, respectively (fig. 2).

. 2.

Relative fitness of the single mutants and populations at 42.2 and 37.0 °C. (A) Relative fitness of the single mutants in competition with the ancestor at 42.2 and 37.0 °C. (B) Relative fitness of the 26 rescue populations in competition with the ancestor at 42.2 and 37.0 °C. Boxplots represent the relative fitness (wr) values of all replicates for each single mutant or population in competition with the ancestor. A wr value near or at 1.0 indicates similar fitness to that of the ancestor; values >1.0 indicate higher fitness than the ancestor, and values <1.0 indicate lower fitness than the ancestor.

In contrast, experiments at the ancestral optimum temperature, 37.0 °C, showed that the single mutants had an average wr disadvantage of 3.5% (one-tailed t-test: P = 0.0002; fig. 2). The rpoB H447L mutant had a particularly low wr value of 0.905, and the hslVU mutations had wr values ranging from 0.933 to 0.993 (table 2). These single mutant results were again nearly identical to results based on population samples, because the 26 rescue populations had a fitness decrease in 3% relative to the ancestor, with 14 of 26 populations having relative fitness values significantly < 1.0 (one-tailed t-test: P < 0.05, fig. 2). Populations with fixed rpoBC mutations had an average wr decrease in 8% (one-tailed t-test: P = 0.0082; n = 7; range: 0.82–0.96), which was significantly lower than all other populations (one-tailed t-test, unequal variance: P = 0.019). Although some hslVU populations had fitness values significantly <1.0 at 37.0 °C (supplementary table S2, Supplementary Material online), hslVU populations had an average wr decrease in 1% relative to the ancestor (n = 16; range: 0.92–1.03), which was not significantly different from 1.0 (one-tailed t-test: P = 0.054).

Gene Expression Differences between Rescue Mutants and REL1206

To attempt to elucidate the molecular mechanisms that lead to evolutionary rescue, we contrasted gene expression between REL1206 and the seven rescue mutants. We gathered replicated RNAseq data for each mutant at two temperatures: 42.2 and 43.0 °C. We used 42.2 °C because it allowed a direct comparison with REL1206 under the same conditions, and 43.0 °C because it is the experimental temperature at which rescue occurred. To assess whether the difference between high stress (42.2 °C) and rescue (43.0 °C) conditions mattered, we performed two analyses. First, we contrasted gene expression between the two temperatures for each mutant. Six of seven mutants had <100 significant DEGs (Padj < 0.001), and the hslU frameshift had 231 DEGs between temperatures (supplementary table S3, Supplementary Material online). Second, we compared the two temperatures for mutants to REL1206 at 42.2 °C. We detected 491 DEGs with the 43.0 °C mutant data and 450 DEGs with the 43.0 °C data, with 75% of DEGs shared between the two analyses (supplementary table S4, Supplementary Material online). Overall, these results suggest some temperature-specific differences between high stress and rescue temperatures. However, comparisons using the mutant data at the two temperatures led to identical trends and qualitative conclusions about functional enrichment and directional changes in gene expression. Hence, for simplicity, we focused on comparisons between REL1206 at 42.2 °C and the single mutants at the rescue temperature —that is 43.0 °C.

At 43.0 °C, single mutations in rpoBC or hslVU exhibited from 58 to 250 upregulated DEGs and between 156 and 369 downregulated DEGs relative to REL1206 at 42.2 °C. We assessed common sets of DEGs among single mutants within specific operons. For example, the five hslVU mutants shared 127 highly downregulated genes (of 294 total) in common, and these exhibited no GO-based enrichment for specific biological processes (fig. 3). Similarly, the two rpoBC mutations shared 250 downregulated genes (fig. 3) that were enriched for transmembrane transporters and catabolic processes (supplementary table S5, Supplementary Material online). Finally, all seven mutants shared 113 downregulated DEGs (fig. 3), a set that could contain genes critical to rescue. GO analyses of this gene set did not reveal a significant enrichment for any biological processes, leaving it difficult to infer which (if any) of these genes contributed to rescue events.

. 3.

The number of upregulated and downregulated genes in the single mutants grown at 43.0 °C relative to the ancestor grown at 42.2 °C. For the Venn diagrams in both A and B, the rpoBC circle represents the number of DEGs shared between the two rpoBC mutants, and the hslVU circle represents genes shared among the five hslVU mutants. (A) The number of downregulated genes relative to the ancestor. (B) The number of upregulated genes relative to the ancestor.

For upregulated genes in the mutants compared with the REL1206 ancestor, we found 185 DEGs shared by rpoBC mutants (of 288 total upregulated genes) and 44 for all five hslVU mutants (of 137 total upregulated genes). Interestingly, 40 of the 44 were also upregulated in the rpoBC mutants (fig. 3). GO analyses on this set revealed enrichment for flagellum assembly and motility (supplementary table S6, Supplementary Material online). More specifically, 16 of these 40 genes were annotated to be directly involved in flagellum regulation, assembly, or motility (Liu and Ochman 2007; Kaundal et al. 2020). Through manual investigation of the remaining 24 genes, we found 11 genes that were involved in membrane transport and 13 genes involved in translational processes, amino acid synthesis, and nucleotide synthesis (supplementary table S7, Supplementary Material online).

Rescue Predominantly Restores Gene Expression

An ongoing question about molecular adaptation is whether it restores physiological and molecular processes from a stressed state back toward the unstressed, wild-type state or whether it instead tends to drive the evolution of novelty (Carroll and Marx 2013). Previous studies have suggested the former, because studies have shown that E. coli adapts to high temperature stress (42.2 °C) by restoring both gene expression (Rodríguez-Verdugo et al. 2016; González-González et al. 2017) and phenotypic characteristics (Hug and Gaut 2015) toward that of the unstressed ancestor.

Following previous studies, we investigated gene expression among mutants, the ancestor at 42.2 °C and the ancestor at 37.0 °C. Similar to those studies, we found a strong negative correlation between ratios that measure the degree of gene expression change in the mutant relative to the two states of the ancestor. For example, the rpoB H447L mutation exhibited a correlation of −0.834 (fig. 4), illustrating that the mutant tended to move gene expression back from the stressed (42.2 °C) state toward the wild-type (37.0 °C) state. Similar negative correlations were obtained with the other seven mutants (fig. 4 and supplementary fig. S1, Supplementary Material online), but the negative correlations were generally stronger for the rpoBC mutations than the hslVU mutations.

. 4.

Direction of gene expression change in single rescue mutants. (A–C) The y axis represents the direction of gene expression change of the three single rescue mutants, (A) rpoB H447L, (B) hslU G60D, and (C) hslV H68P, when grown at 43.0 °C compared with the ancestor grown at 42.2 °C. These three mutations were chosen as illustrative, with the remaining four mutations shown in supplementary figure 2, Supplementary Material online. The x axis represents the ancestral changes in gene expression when it is grown at 42.2 °C compared with when it is grown at 37.0 °C. The black line represents the linear regression fitted to the data in each graph. (D) Comparison of the number of genes in each category of gene expression change for the rescue mutants and for the rpoB and rho mutants studied in Rodríguez-Verdugo et al. (2016) and González-González et al. (2017). The rpoB and rho mutations were adaptive to high but nonlethal temperatures.

We also counted the number of genes that fell into one of four expression categories: restored, unrestored, reinforced, or novel (see Materials and Methods section). The predominant category was restored, which suggests that the mutants shifted gene expression back toward the unstressed state (table 3 and fig. 4). Five of seven mutants had >50% of their genes in this category, whereas the hslU frameshift mutant and the hslU G60D mutants had 38% and 40% of their genes restored, respectively. The next highest category was unrestored, which represented 23–53% of the genes in the rescue mutants. Perhaps the most striking aspect of our analyses was that each rescue mutant had >100 genes that exhibited novel expression patterns. This category had far fewer genes in previous studies; for example, high temperature adaptive mutations in rpoB and rho caused <60 genes to have novel expression patterns (Rodríguez-Verdugo et al. 2016; González-González et al. 2017). Given this apparent difference, we compared the average number of novel genes across all seven mutations to previous studies that used the same methods (Rodríguez-Verdugo et al. 2016; González-González et al. 2017). We found that rescue mutations had a significantly different proportion of genes in the four categories to previously studied rpoB and rho mutations (P < 2.2 × 10−16, contingency test; fig. 4). A total of 43 genes displayed novel expression patterns across all seven rescue mutants; according to GO analyses, these were enriched for transmembrane transporters for carbohydrates, such as glucose and mannose (supplementary table S8, Supplementary Material online).

Table 3

Number of Genes in Each Category of Gene Expression Change

	Restored	Reinforced	Unrestored	Novel	Total
hslU frameshift	669	8	923	135	1,735
hslU G60D	719	9	872	183	1,783
hslU L163R	960	4	636	121	1,721
hslV H68P	947	7	646	154	1,754
hslV frameshift	977	4	619	118	1,718
rpoB H447L	1,144	7	449	305	1,905
rpoC W1020G	1,077	3	520	159	1,759

Discussion

We have performed experiments to characterize the genetic mutations that contribute to the rescue of E. coli populations from an otherwise lethal temperature of 43.0 °C. Overall, we have found that rescue is infrequent, because it occurred for only 8.8% (26 of 296) of our experimental populations. Although not all of our populations were independent (fig. 1 and Materials and Methods), the observed rescue frequency is similar to that of Mongold et al. (1999), who found that 10% of E. coli populations recovered from a 44.0 °C treatment (Mongold et al. 1999). One difference is that they observed rescue only in populations that were preadapted to thermal stress, whereas our populations were not preadapted. Nonetheless, the two studies are consistent in showing that evolutionary rescue is infrequent but, somewhat paradoxically, frequent enough to be a potent source of evolutionary innovation (Bell 2017).

We have studied thermal stress because it has complex effects on a wide variety of physiological functions, implying that rescue adaptations could be genetically diverse. Among the 26 populations that exhibited U-shaped rescue dynamics (fig. 1), we have identified 29 distinct fixed mutations (table 1). There were clear patterns among these mutations, because some mutations were found in parallel across presumably independent experiments. For example, the same point mutation in ECB_00530 was found in 3 separate weeks, as was a nonsynonymous mutation in rpoC (table 1). Other common locations of fixed mutations included the hslVU heat shock protease operon, the rpoBC RNAP operon, and an intergenic region between clpA and serW, which encode a component of a protease system similar to that encoded by hslVU (Kwon et al. 2004) and a serine-bearing tRNA.

Mutations in some of our genes have been identified in previous experiments of E. coli temperature adaptation under nonlethal conditions (Tenaillon et al. 2012; Deatherage et al. 2017). For example, mutations in mrdA and rpoBC have been identified in numerous evolution experiments, both in low-nutrient conditions and under temperature adaptation (Conrad et al. 2010; Tenaillon et al. 2012; Long et al. 2015; Deatherage et al. 2017). Similarly, multiple mutations in hslVU were identified in an experiment that evolved REL1206 for 2,000 generations at several temperature regimes, including 37.0 and 42.0 °C (Deatherage et al. 2017). However, the hslVU mutations were primarily observed at 37.0 °C, not at the stressful temperature, and the 37.0 °C mutations did not obviously inhibit function via frameshifts or premature stop codons. Interestingly, mutations in hslVU were not found commonly during evolution at 42.2 °C by Tenaillon et al. (2012). Of their 115 lines and >1,000 mutations, only one line carried a nonsynonymous mutation in hslU. These observations suggest that slight differences in conditions (i.e. from stressful to lethal temperature) may have large effects on the set of potentially adaptive mutations.

Single Mutations Drive Rescue

Of our 26 rescued populations, 10 had only a single fixed variant, suggesting they were drivers of evolutionary rescue. To explore the fitness effects and the potential mechanistic basis of these potential drivers, we isolated single mutant genotypes for seven mutations within the rpoBC and hslVU operons (table 2). We isolated two nonsynonymous mutations in the former, which encodes the β and β′ subunits of RNAP. Mutations in RNAP must maintain enzyme function due to its central role in transcription, but it is also known that nonsynonymous mutations in RNAP can have numerous effects on cellular properties like fitness, growth rate, and patterns of gene expression (Herring et al. 2006; Rodríguez-Verdugo et al. 2014, 2016; Carroll et al. 2015).

We also isolated three nonsynonymous mutations and two frameshifts within the hslVU operon (table 1), which encodes two heat shock proteins that form an ATP-dependent protease complex. In hslU, one nonsynonymous mutation (G60D) is in the N-terminal domain that has ATPase activity, the other (L163R) is in the I-intermediate domain that recognizes substrates, and the frameshift is in the C-terminal domain that interacts with the hslV protein product (Bochtler et al. 2000; Lien et al. 2009). Previous studies have detailed the effects of nonsynonymous mutations in the N- and C-terminal domains of hslU and concluded that most mutations cause the loss of ATP hydrolyzing ability and protease activity (Shin et al. 1996; Bochtler et al. 2000). Similarly, research has shown that nonsynonymous mutations throughout the sequence of hslV causes reduced protease activity (Yoo et al. 1997, 1998). Taken together, these studies suggest that the hslVU mutations in our study likely reduce or completely knockout their heat shock protease activity.

We assayed the relative fitness of each of seven mutations to confirm that they can drive rescue dynamics. Six of the seven have wr values significantly greater than 1.0 at high temperature, with the last borderline significant (P = 0.073). The wr estimates range from a ∼6% fitness increase for hslU G60D to 30% or higher fitness increase for three of the seven mutations (the hslU frameshift, hslU L163R, and rpoB H447L; table 2 and fig. 3). This range of wr values is not particularly unexpected, even for mutations within the same gene (Barrick et al. 2010; Conrad et al. 2010; LaCroix et al. 2015). For example, Rodríguez-Verdugo et al. (2014) assessed the fitness of four adaptive mutations in two different codons of rpoB, and their fitness benefits varied from 17% to 37%. Similarly, González-González et al. (2017) found that adaptive mutations within the rho gene varied in fitness increases from 8% to 26%. Together, the wr estimates of our seven mutations, coupled with the fact that each was the lone fixed mutation in at least one rescue population, clearly establish that each mutation is sufficient for rescue.

Population Dynamics of Rescue Mutations

There were, however, at least two interesting differences in the observed patterns of hslVU and rpoBC mutations. First, fixed mutations in these two operons were not found together in the same population, which is statistically improbable given their respective frequencies across populations (P < 0.02). The lack of co-occurring rpoBC and hslVU mutations suggests that clonal competition canalizes the initial adaptive response. To the extent that wr values at 42.2 °C reflect fitnesses at 43.0 °C, the relative fitness assays suggest that the two rpoBC mutations would outcompete at least three of the hslVU mutations when both are present (fig. 2). Second, patterns of parallelism differed between hslVU and rpoBC mutations. Specific hslVU mutations were found across different weeks and among nonindependent populations within weeks (fig. 1). In contrast, rpoBC mutations were identified across weeks but typically in only a single population.

What might drive these apparently different patterns? To address this question, we first recognize that evolutionary rescue can act on standing genetic variation prior to the introduction of the lethal stressor (Bell 2013, 2017). Indeed, our experiments were unlikely to be severely mutation-limited (Lang et al. 2013) because the experiment for each week began in an overnight DM25 culture at 37.0 °C (fig. 1). We nonetheless believe two factors may have contributed to different patterns for hslVU and rpoBC rescue mutations. First, there is strong constraint on function for rpoBC mutations, whereas hslVU knockouts are adaptive at lethal temperatures. Hence, we suspect that there are more potential rescue mutations in hslVU due to fewer functional constraints. This difference may alone explain why fixed hslVU mutations were more numerous than rpoBC mutations (table 1), but it does not fully explain why hslVU mutations tended to be more common across populations in 1 week.

Second, we suspect that trade-offs contribute to the pattern across populations in 1 week. We posit that fitness costs at 37.0 °C affect the frequency of individual mutations in the overnight culture, which in turn affects the probability of a mutation being sampled into multiple populations in any given week (fig. 1). We thus assessed wr for the single mutants at 37.0 °C, the temperature of the initial batch cultures, for trade-offs. We found that fitness costs are not a universal feature of the rescue mutations, at least within the power of our experiments to detect such differences. Only two of the seven mutations exhibited significant fitness deficits: rpoB H447L and the hslU frameshift mutation (table 2). These observations contribute to a growing consensus that trade-offs, and specifically antagonistic pleiotropy, are common but not universal (Cooper and Lenski 2000; MacLean et al. 2004; Rodríguez-Verdugo et al. 2014; Deatherage et al. 2017). The pattern of trade-offs does not fully support our model, because the hslU frameshift was common across populations in week 3 despite its low fitness at 37.0 °C (table 1). Nonetheless, wr values based on populations consistently show that the populations with fixed rpoBC mutations do exhibit trade-offs (fig. 2). We thus continue to suspect that trade-offs play a large role in the population dynamics of our experiment, because mutations with trade-offs are at low(er) frequency in the 37.0 °C overnight culture and thus less likely to be sampled into multiple populations.

We add two additional points about the dynamics of rescue in our experiment. First, the fact that we uncover clear and repeatable patterns of mutations in only a small subset of genes and operons suggests that the universe of potential rescue mutations is small, especially given that the experiment should not have been severely mutation-limited. One potential explanation for parallel mutations across weeks—such as the clpA/serW mutation and several others (table 1)—is that the experiment selected for low-frequency mutations that were present in the common, frozen ancestral stock (fig. 1). However, this possibility does not contradict (but rather reinforces) the conjecture that there are only a few major mutational targets for adaptive rescue. In this context, it is interesting to muse whether there is in fact a very small universe of mutations that are capable of rescue or whether there is a large universe of such mutations but most do not establish in our populations because they have severe fitness costs at 37.0 °C. We cannot yet distinguish between these two alternatives, but knowing the prevalence of trade-offs is important to furthering our understanding about the dynamics of evolutionary rescue. Second, we need to mention an important distinction between our E. coli experiment and evolutionary rescue in natural populations, particularly in size-limited, nonbacterial populations subjected to stressors like climate change. The yeast experiment under lethal salt conditions is illustrative, because it showed that rescue occurs less frequently in populations of small size (Bell and Gonzalez 2009). It thus seems likely that the frequency of rescue in most plant and animal populations, which tend to be relatively small, will be far less than the 8–10% estimated for E. coli under lethal temperature stress.

Insights into Rescue Mechanisms and Evolutionary Direction

We have established that the seven mutations are capable of rescue. However, an overarching question is about their function—that is, what do they do and how do they drive rescue events? To begin to address this question, we generated gene expression data for the seven single rescue mutants. For each clone we measured gene expression at the exponential phase of growth at two temperatures (42.2 and 43.0 °C) to compare with the REL1206 ancestor at 42.2 and 37.0 °C. Our goals were: 1) to find sets of DEGs in common across the entire set of mutants, in the hope that they yield clues to mechanisms and 2) to characterize the overall direction of gene expression changes with respect to the stressed and unstressed state of the ancestor.

Our first goal was formulated under the hypothesis that rescue mutants may affect common pathways that lead to rescue. We proffer this hypothesis knowing that the differences in wr among mutants reflects the fact that many DEGs vary among them and also that different mutants may have utilized different pathways to achieve rescue. Nonetheless, we first focused on shared DEGs between the two rpoBC mutants. They share a set of 435 DEGs relative to the ancestor (fig. 3), based on a conservative measure of differential expression that includes both adjusted P values <0.001 and 2-fold differences in log2 expression (see Materials and Methods section). The high number of DEGs is not surprising, because previous work has shown that single mutations in RNAP can alter the gene expression of ∼1,000 or more genes (Conrad et al. 2010; Carroll et al. 2015; Rodríguez-Verdugo et al. 2016). The number of common DEGs was much lower for hslVU mutations, at 171 total; the union between the two sets yielded 40 up- and 113 downregulated genes (fig. 3).

The set of 113 downregulated genes provided no clear patterns with regard to function, based on GO analyses and manual investigation of gene annotations. However, the 40 upregulated genes yielded two notable observations. First, 24 of the 40 were involved in transport, ribosomal assembly, and amino acid and nucleotide pathways. This suggests that changes in gene expression were supporting translation processes, perhaps to enhance efficiency at high temperature. Surprisingly, GO analyses revealed that this 40 gene set is also enriched for genes involved in flagellar assembly and motility, a process that is known to be energetically costly (Soutourina and Bertin 2003). We manually verified that at least 16 of the 40 genes have functions associated with flagella. Given this observation, one must ask about the potential adaptive benefit of flagellar production and activity. It is hard to answer this question directly, but it has been shown that E. coli produce flagella in low nutrient environments (Shi et al. 1992; Sim et al. 2017). Hence, one hypothesis is that increased expression of flagella leads to enhanced motility and nutrient acquisition. It seems doubtful, however, that enhanced motility is an advantage in our well-mixed system. Interestingly, previous work has shown that first-step adaptive mutations also increased flagellar gene expression, only to be attenuated by subsequent compensatory mutations (Rodríguez-Verdugo et al. 2016). These observations suggest that flagellar production may be a disadvantageous by-product of other major and adaptive shifts in physiological processes. Altogether, then, it is not clear which—if any—of the common DEGs contribute to the rescue phenotype. We suspect, however, that shifts in the expression of translation-related genes are more critical for adaptation than the upregulation of flagella-related genes. In the future, proteomic analyses may provide further insights into changes introduced by rescue mutations and evolutionary mechanisms.

It is worth briefly considering the potential effects of hslVU mutations separately, because the frameshift mutations lead to the somewhat paradoxical conclusion that knockouts of heat shock-related proteins are beneficial at lethal temperatures. In this context, it is helpful to know that the hslVU protein degrades the σ32 factor. In E. coli, σ32 typically exists as an RNA with secondary structure that unfolds under high temperature to allow translation, and then the σ32 protein regulates the transcription of genes needed to carry out the heat shock response (Roncarati and Scarlato 2017). By degrading σ32, the hslVU protein indirectly inhibits the production of other proteins in the heat shock cascade (Kanemori et al. 1997). We propose that the hslVU rescue mutations decrease or knockout the protease function, thereby facilitating an uninhibited heat shock response through σ32. One interesting fact is that σ32 regulates a series of genes—like dnaK, dnaJ, and grpE (Nonaka et al. 2006)—that are also required for the control of flagellar synthesis through another sigma factor (σ28) (Shi et al. 1992), suggesting a mechanistic link between hslVU and flagellar genes. These observations support our hypothesis that loss of hslVU activity enhances some aspects of the heat shock response and strengthens the possibility that the upregulation of flagellar genes is not a direct feature of adaptation.

In a second goal, we used expression data to characterize the directional change for all genes (Carroll and Marx 2013). The question is whether our rescue mutations produce patterns similar to previously studied mutations that contributed to temperature adaptation (Rodríguez-Verdugo et al. 2016; González-González et al. 2017). The previous mutations predominantly moved gene expression from a stressed physiological condition toward the unstressed, wild-type condition. Our rescue mutations are similar, because they also predominantly shift expression toward restoring the unstressed state (fig. 4). However, there is an important difference: the rescue mutations yielded significantly more genes with novel gene expression patterns (fig. 4), particularly for transmembrane transporters. This result is not dependent on the 43.0 °C conditions, because it is also evident at the slightly modified temperature of 42.2 °C, under conditions identical to the previous experiments (supplementary figs. S2 and S3 and table S9, Supplementary Material online for 42.2 °C results).

These data suggest that there could be a qualitative difference between rescue mutations and mutations that contribute to adaptation under nonlethal conditions. This conclusion is of course subject to caveats. For example, although we followed procedures identical to previous studies, the differences may nonetheless reflect experiment-specific effects, such as batch effects in RNAseq data (although we explicitly examined such effects; see Materials and Methods section), rather than differences in the dynamics of adaptation. We also assayed only a single point in the growth phase, which may not represent the crucial point in the cell cycle for rescue adaptations. As a consequence, the novel expression patterns we have observed may represent noise rather than changes that contribute to (or are necessary for) adaptation. However, it is possible that phenotypic novelty is an important feature of evolutionary rescue; that is, when challenged with a lethal stressor, it is not sufficient to move toward restoring expression toward the unstressed state. Ours is only a first observation, but it opens an interesting question for future research: do mutations that drive evolutionary rescue differ qualitatively from adaptations in nonlethal environments?

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Data Availability

All high-throughput sequence data generated in this study have been submitted to the NCBI short-read archive at https://www.ncbi.nlm.nih.gov/sra and can be accessed with project numbers PRJNA326455 and PRJNA640188.

Click here for additional data file.