Phenotypic and Genotypic Adaptation of Escherichia coli to Thermal Stress is Contingent on Genetic Background

Abstract Evolution can be contingent on history, but we do not yet have a clear understanding of the processes and dynamics that govern contingency. Here, we performed the second phase of a two-phase evolution experiment to investigate features of contingency. The first phase of the experiment was based on Escherichia coli clones that had evolved at the stressful temperature of 42.2 °C. The Phase 1 lines generally evolved through two adaptive pathways: mutations of rpoB, which encodes the beta subunit of RNA polymerase, or through rho, a transcriptional terminator. We hypothesized that epistatic interactions within the two pathways constrained their future adaptative potential, thus affecting patterns of historical contingency. Using ten different E. coli Founders representing both adaptive pathways, we performed a second phase of evolution at 19.0 °C to investigate how prior genetic divergence or adaptive pathway (rpoB vs. rho) affects evolutionary outcomes. We found that phenotype, as measured by relative fitness, was contingent on founder genotypes and pathways. This finding extended to genotypes, because E. coli from different Phase 1 histories evolved by adaptive mutations in distinct sets of genes. Our results suggest that evolution depends critically on genetic history, likely due to idiosyncratic epistatic interactions within and between evolutionary modules.


Introduction
famously argued that historical contingency is a defining feature of evolution. He asserted that natural selection is deterministic in the absence of historical contingency, and hence it can, in theory, converge on an optimal solution to any challenge. In contrast, evolution may be contingent on the idiosyncratic nature of historical events, leading to unpredictable and perhaps unrepeatable evolutionary outcomes. Historical contingency can also reflect "causal dependance", the notion that the outcome of evolution depends on the preceding state, as opposed to different outcomes resulting from the same state (Beatty and Carrera 2011).
It is obvious that evolution must be historically contingent to some degree, but questions remain about the patterns and mechanisms of contingency. Exploring the interplay of contingency and determinism, and the genetic effects that drive their dynamics, is crucial for understanding the evolutionary process. Moreover, understanding this interplay is helpful for predicting outcomes to applied problems-for example, projecting which species will survive climate change, identifying evolutionary constraints on evolving resistance to chemotherapeutic agents, and forecasting pathogen variation, mutation, and epidemiology (Bay et al. 2017;Vlachostergios and Faltas 2018;Leray et al. 2021).
Historical contingency has been studied in both natural populations and in the laboratory. In the field, experiments have inferred deterministic outcomes from convergent evolutionary events (Losos 2011). For example, one experiment tested brown anole lizard populations that were subjected to living on narrow perches and found that all lizard populations evolved shorter limbs (Kolbe et al. 2012). Similarly, male guppies across different populations evolved shorter life histories in the absence of predators (Reznick and Bryga 1987). These pervasive convergent outcomes have been used to argue that natural selection is predictable and that historical contingency does not play a major role in evolution (Losos 2010;McGhee 2016). This interpretation is debated; however, because natural populations may contain standing genetic variation that increases the probability of parallel responses (Blount et al. 2018).
In the laboratory, controlled experiments have investigated contingency by evolving populations from a single ancestral genotype (Blount et al. 2018). For example, in the long-term evolution experiment (LTEE), 12 replicated Escherichia coli populations converged phenotypically, as evidenced by increased fitness, faster growth, and larger cells (Bennett and Lenski 1993;Wiser et al. 2013), suggesting a lack of contingency based on mutations that arose during the experiment (Blount et al. 2018). However, Wiser et al. (2013) documented significant variation in the fitness trajectories among LTEE lines, suggesting evolutionary trajectories are influenced by history despite the fact that overall phenotypic outcomes appear to converge. Additionally, one rare adaptation in the LTEE, aerobic growth on citrate (Blount et al. 2008), was contingent on previously occurring, "potentiating" mutations, illustrating causal dependance. Here, antecedent genetic changes were essential for phenotypic innovation and were co-opted for future innovative leaps (Quandt et al. 2015). Further work showed that aerobic citrate utilization evolved via a genetic trajectory that maintained the potential for evolving the phenotype (Leon et al. 2018), again highlighting the influence of history on evolution.
Another example from the LTEE is the evolution of antibiotic resistance, because both the level of resistance and the complement of resistance mutations varied as a function of the starting genetic background of the LTEE line (Card et al. 2021). These studies show that the evolution of phenotypes can be, but are not always, historically contingent on the genetic background. Nonetheless, important questions remain (Blount et al. 2018). For example, what circumstances favor contingent versus deterministic evolutionary outcomes? How does prior genetic divergence between populations affect the likelihood of parallel evolutionary responses? And to what extent does epistasis shape evolutionary outcomes?
One approach to study these questions is two-phase (or "historical difference"; Blount et al. 2018) evolution experiments ( fig. 1). In this experimental set-up, Phase 1 consists of initially identical populations that are evolved in the same environment, leading to potential genetic and phenotypic divergence among the populations. These populations are then subjected to a second phase (Phase 2) by allowing them to evolve in a new environment. The question is whether Phase 2 evolution differs across populations as a function of Phase 1 history. That is, is evolution in Phase 2 constrained by (or contingent upon) phenotypic and genotypic variation that was generated in Phase 1?
Thus far, the two-phase approach has uncovered nuanced insights into contingency. For example, Plucain et al. (2016) evolved 16 E. coli populations in four different chemical environments for 1,000 generations (Phase 1) before propagating them for another 1,000 generations in a single, new environment (Phase 2). The evidence for contingency was mixed: they found some evidence for contingency on phenotypic evolution, because both the growth rate and fitness of Phase 2 populations varied according to their Phase 1 history. However, they found no evidence for contingency at the molecular level; the mutations that arose during Phase 2 evolution did not vary as a function of the genetic backgrounds generated during Phase 1 (Plucain et al. 2016). Another two-phase evolution experiment used yeast and estimated that ∼50% of the variance in fitness across Phase 2 populations was attributable to Phase 1 history, mostly because Phase 1 populations with low fitness evolved more rapidly in Phase 2 (Kryazhimskiy et al. 2014). However, like the E. coli study, this study also found no evidence to suggest that Phase 2 genotypic changes were influenced by Phase 1 history.
Based on their two-phase experiment, Kryazhimskiy et al. (2014) proposed the "global epistasis hypothesis" a form of diminishing returns epistasis. Diminishing returns epistasis implies that adaptive mutations have larger selective effects in relatively unfit genotypes (Griffing 1950;Jerison and Desai 2015). The global epistasis hypothesis further posits that a mutation's effect depends solely on the fitness of the genetic background and not on the background genotype (Wei and Zhang 2019). In this framework, evolutionary trajectories are predictable based on fitness information alone. In contrast, recent work has suggested that epistasis may be idiosyncratic, in that the direction and/or magnitude of epistatic interactions depends on specific genotypes (Wei and Zhang 2019;Bakerlee et al. 2022), potentially making evolutionary outcomes less predictable. Idiosyncratic epistasis may even be modular, because it is a property of interactions among genes that contribute to specific functions (Tenaillon et al. 2012) or that vary by environment (Wei and Zhang 2019).
To examine these ideas further, we perform a two-phase evolution experiment in E. coli that uses extreme temperatures as the selective environments. We focus on extreme temperatures for four reasons. First, temperature is a fundamental environmental property that affects physiological traits and often defines species' distributions; hence, it often requires a complex evolutionary response (Somero 1978;Cooper et al. 2001). Second, temperature adaptation often, but not always, leads to trade-offs at other temperatures (Rodríguez-Verdugo et al. 2014), suggesting that contingency could be important in this system. Third, temperature is a topic with rich historical precedent in the experimental evolution literature. For example, Bennett and Lenski (1993) evolved E. coli at 20 °C, near the lower edge of the temperature niche, after first adapting them to the upper end of the temperature niche (42.2 °C). They found no convincing evidence of contingency at the phenotypic level, but their work was based on a relatively small number of samples (n = 12) and lacked genetic information.
The fourth reason that we focus on temperature is that we can take advantage of a previous large-scale experiment. Tenaillon et al. (2012) evolved 115 lines from a single E. coli founder strain (REL1206) at 42.2 °C. After 2,000 generations of evolution, they evaluated a single clone from each population for fitness gains and sequence changes. The sequence changes revealed that adaptation often occurred through two distinct adaptive pathways defined by mutations in either the RNA polymerase subunit beta (rpoB) gene or the transcriptional terminator (rho) gene. Mutations in both of these genes occurred statistically less often than expected by chance, suggesting negative epistatic interactions. Intriguingly, the two pathways were each positively associated with additional distinct sets of mutations. For example, rpoB clones tended to have mutations in rod, ILV, and RSS genes, but mutations in these genes were rare in the rho lines. The complex landscape of both negative (e.g., rpoB vs. rho) and positive (e.g., rpoB with rod and ILV) associations suggests that evolutionary changes are partially dependent on the genetic background and that the two pathways may represent discrete evolutionary modules.
Here, we hypothesize that these associations constrain future adaptation and thus affect patterns of historical MBE contingency. To test this hypothesis, we utilize evolved clones from Tenaillon et al. (2012) to represent Phase 1 of a two-phase experiment. After choosing a set of clones representing rpoB and rho genotypes, we evolve them at the lower extreme of E. coli's temperature niche (19.0 °C) and then measure phenotypic and genotypic differences among evolved populations. In doing so, we address three sets of questions: First, do the rpoB and rho lines differ in their response to selection at 19.0 °C, as measured by their fitness response? We are particularly interested in testing one of the predictions of the global epistasis model, which is that populations with lower fitness should exhibit larger fitness gains. Second, is there evidence to suggest that the evolution of rpoB and rho lines differ in their genotypic patterns of change? That is, do the mutations that appeared during Phase 1 evolution shape the set of adaptive mutations that accumulate in Phase 2? Finally, what do our results imply about the evolutionary process, particularly whether epistasis is global or more idiosyncratic and modular?

Selecting Phase 2 Founders Representing two Adaptive Pathways
Phase 1 consisted of 114 lines evolved at 42.2 °C (Tenaillon et al. 2012). After 2,000 generations of evolution, single clones from these lines experienced fitness gains of ∼42%, on average, relative to the single founding ancestor, which we call the Phase 1 Ancestor ( fig. 1). To perform our Phase 2 experiment, we selected five rpoB clones and five rho clones as the founding genotypes for evolution at 19.0 °C. We refer to these ten clones as the Phase 2 Founders ( fig. 1) and label each by its founding line and by its rho or rpoB genotype (e.g., rho_A43T) (table 1). It is important to recognize, however, that the Phase 2 Founders had mutations in additional genes (i.e., not just rpoB and rho) relative to the REL1206 Phase 1 Ancestor (Tenaillon et al. 2012; supplementary table S1, Supplementary Material online). For example, six of the Phase 2 Founders had a mutation in the cls gene, and four had a mutation in ybaL, although mutations in neither gene were strictly associated with either of the two adaptive pathways.
The  1. A schematic of the two-phase evolution experiment. The first phase of evolution was described in Tenaillon et al. (2012). The second phase used a subset of evolved clones from Phase 1 and represented three pathways: clones with rpoB mutations (five genotypes), clones with rho mutations (five genotypes), and the Phase 1 Ancestor.

MBE
and 0.990 for rho (t-test, P = 0.095). We note, however, that the w r variance was higher among the rpoB founders (Var(w r ) = 0.0033) compared with the rho founders (Var(w r ) = 0.00075). Fourth, we selected clones with a single mutation in either rpoB or rho but not in both genes. Finally, we chose a sample of distinct rpoB and rho mutations, because Phase 1 mutations occurred across different codons and caused different amino acid replacements (table 1).
Once chosen, Phase 2 Founders were propagated at 19.0 °C for 1,000 generations, with six replicate populations per founder, under conditions identical to the Phase 1 experiment except for temperature (19.0 °C vs. 42.2 °C). We also evolved 12 replicates of the Phase 1 Ancestor as a control ( fig. 1), making a total of 72 (=6 × 10 + 12) populations in the Phase 2 experiment. Among the 72 populations, 7 went extinct, including 1 of the 12 descended from the Phase 1 Ancestor, 4 of 6 populations descended from Phase 2 Founder Line 3 (rpoB I966S), and 2 populations descended from Phase 2 Founder Line 142 (rpoB I572L) (table 1). The following analyses were therefore performed on the set of 65 surviving populations.
Relative Fitness at 19.0 °C Varies Significantly Among Pathways and Founder Genotypes Previous two-phase experiments have shown that fitness can be affected by historical contingency (Kryazhimskiy et al. 2014;Plucain et al. 2016). To test for such contingencies, we first measured w r of the Phase 2 populations against the Phase 1 Ancestor, based on three technical replicates per population. From a total of ∼200 competition experiments, we estimated that the fitness of evolved populations increased by 3.6%, on average, at 19.0 °C (P < 0.01, Wilcoxon test; fig. 2). On average, lines descended from rpoB backgrounds had 1.0% higher fitness than the Phase 1 Ancestor; the rho lines had higher fitness by 6.4%; and the control lines increased by 3.4% (table 2). To test whether these differences were significant, we applied Analysis of Variance (ANOVA) that partitioned by pathway (rho vs. rpoB) and were nested by Phase 2 Founder genotypes (table 2). The pathway effect was significant (P < 0.003) and explained 10.0% of the w r variance, but the Phase 2 Founder genotype explained an even higher proportion of the variance (23.0%; P = 0.11). We also applied Analysis of Covariance (ANCOVA) using the fitness of Phase 2 Founders as covariates, and the difference between pathways remained significant (P = 0.029, ANCOVA). Finally, we note that w r values were not normally distributed (Shapiro-Wilk, P < 0.01), despite the large sample size and even after routine normality transformations, thus violating assumptions used in ANOVA and ANCOVA. However, we also repeated the analysis using nonparametric tests and obtained similar results (Kruskal-Wallis; effect of pathway: P = 0.0002, η² = 0.0784; effect of Phase 2 Founder: P = 0.0002; η² = 0.153). Overall, these observations indicate that the fitness response depended on the pathway but also on the Phase 2 Founder genotypes within pathways.
Previous work has shown that the rate of change of Phase 2 populations can vary as a function of the fitness of Phase 2 Founders-that is, less fit Founders led to generally larger leaps in fitness during Phase 2 evolution (e.g., Kryazhimskiy et al. 2014). To assess this potential effect, we first measured the fitness of Phase 2 populations by competing them against their respective Phase 2 Founders at 19.0 °C, constituting another ∼160 w r competition assays (table 2). On average, the Phase 2 populations had w r = 1.08, thus reflecting, on average, a significant 8% fitness advantage at the end of the experiment (P < 2.20 × 10 −16 , one-sample t-test) compared with their Founders. Populations descended from rpoB genotypes experienced a 9% fitness advantage on average (P = 1.41 × 10 −15 , onesample t-test), those descended from rho founders had a

MBE
7% fitness advantage (P = 1.32 × 10 −10 , one-sample t-test), and the difference was not significant (P = 0.063, unpaired t-test; fig. 3A). Eight of the 10 sets of Phase 2 populations had significant fitness advantages relative to their Founders (table 2 and fig. 3B), and the effect of the Phase 2 Founder genotype on fitness was again significant (η² = 0.23, P = 1.49 × 10 −6 , Kruskal-Wallis test). We then plotted w r for all Phase 2 Founders against the difference in w r between the Phase 2 Founder and each of its evolved populations ( fig. 3C). (These w r values were all

FIG. 2.
Relative fitness (w r ) of the Phase 2 evolved populations at 19.0 °C based on competition assays to REL1207 (an Ara + variant of the Phase 1 Ancestor; see Materials and methods section). On the x-axis, Control represents the Phase 1 Ancestor competed against itself and shows that the Phase 1 Ancestor (REL1206) and the REL1207 variant have similar fitnesses, as expected. The boxplot for Phase 1 Ancestor reports w r for the control populations evolved from REL1206, and the remaining boxplots show w r values for different Phase 2 founders. The dots in each boxplot show each w r measurement across replicated populations, with each population measured at least three times. The horizontal line within the boxplot is the median w r value, with the box representing the upper and lower quartile and whiskers are calculated using the interquartile range. MBE measured relative to the Phase 1 Ancestor.) If less fit Founders lead to generally larger fitness gains, we expected a negative slope. The combined set of rpoB and rho populations followed this prediction because lower-fitness founders had slightly bigger shifts in fitness during Phase 2 evolution (slope = −0.87, P = 0.19, fig. 3C). Individually, the trend was especially evident among the rpoB populations (rpoB slope = −1.32, P = 0.014; fig. 3C), but it did not hold for rho Founders and descendant populations. In fact, the slope based on rho populations was positive, although not significantly so (rho slope = 1.01, P = 0.14; fig.  3C), such that the rho and rpoB slopes were statistically different (P < 0.001). Although the cause(s) of these different patterns between pathways was not clear, it suggests, at a minimum, that the rate of fitness change during Phase 2 evolution was not a simple function of the fitness of Phase 2 Founders.
High-Temperature Trade-offs are Contingent on the founders' Adaptive History and Genotype Previous research has demonstrated significant differences in trade-off dynamics between high-temperature adapted genotypes. For example, nearly half of the 42.2 °C adapted lines from Tenaillon et al. (2012) were less fit than the Phase 1 Ancestor at lower temperatures (37 °C and 20 °C), slightly more than half exhibited no obvious tradeoff, and a surprising few were actually fitter than the ancestor at low temperatures (Rodríguez-Verdugo et al. 2014). These trade-off dynamics imply the possibility of genotype-specific contingencies; we thus investigated trade-off dynamics at 42.2 °C for Phase 2 populations relative to their Phase 2 Founders. As expected, the Phase 2 evolved populations generally had lower fitness (average w r = 0.89, P = 1.14 × 10 −14 , Wilcoxon test) than their Founders at 42.2 °C. Of the ten rho and rpoB Phase 2 groups, seven of ten had significantly lower w r at 42.2 °C compared with their Phase 2 Founder (table 2). We investigated whether these patterns mapped to either pathways or starting genotypes; the difference in average w r between rho and rpoB populations was not statistically significant (P = 0.38, Wilcoxon test; fig. 3D). It was nonetheless notable that rho lines experienced a fitness decline (relative to their Founder) of 8.5% whereas rpoB declined by 15% on average, suggesting some difference in trade-off dynamics between pathways. Moreover, the starting genotype had a significantly large effect on the values of w r at 42.2 °C (η² = 0.34, P = 1.36 × 10 −9 ,

MBE
Kruskal-Wallis test; fig. 3E). For example, lines descended from the rpoB G556S and rho V206A Founders had extremely low fitness at 42.2 °C, at w r = 0.52 (P < 0.01, Wilcoxon test) and w r = 0.80 (P < 0.01, Wilcoxon test). The take-home point is that trade-offs differed as a consequence of Phase 2 Founder genotypes.
Mutations that Arose During Phase 2 Evolution are Contingent on Genetic History To examine contingency at the level of individual mutations, we sequenced the DNA of all 65 Phase 2 and control populations. After filtering the sequencing data and calling genomic variants, we identified 1,387 point mutations and short indels (<50 bases in length) that arose during the Phase 2 experiment and had population frequencies >5% (supplementary fig. S2, Supplementary Material online). Almost half (45%) of the 1,387 mutations were present at a frequency of <10%, but 119 were "fixed" (as defined by a frequency >85% in a single evolved population). Fixed mutations were present in 51 of the 65 sequenced populations, but there was no discernible pattern across the 14 lines that lacked fixed mutations by pathway (rho, rpoB, and control Founders). Overall, the largest proportion of mutations, 54.4% (742/1387), occurred in intergenic regions ( fig. 4A), and 95.8% of these were point mutations. Within genes, most (89.8% or 371/413) point mutations were nonsynonymous.
We used the sequencing data to confirm that our populations were not cross-contaminated during the 1,000 generation experiment by first assessing whether mutations from the Phase 2 Founder were fixed in the evolved populations, as expected. It was true in every case. We then built phylogenies based on all of the sites in the Phase 2 populations that differed from the Phase 1 Ancestor. The data confirmed expected phylogenetic relationships based on the experimental design (supplementary fig. S3, Supplementary Material online) and thus yielded no evidence of contamination.
Given a lack of obvious evidence for contamination, we first asked whether the patterns and numbers of mutations differed significantly among pathways. In terms of mutations counts, we identified an average of 22 mutations in populations descended from rpoB lines and an average of 20 mutations in populations descended from rho lines that were at a frequency of 5% or higher in the population, a difference that was not significantly different (P = 0.17, unpaired t-test). There was also no difference between pathways in the number of fixed mutations per population (P = 0.14, unpaired t-test). We also contrasted the proportions of mutational variant types (intergenic, frameshift, nonsynonymous, and synonymous mutations and large deletions >50 bps) between rho and rpoB pathways, again finding no difference for the complete set of mutations (P = 0.78, contingency test; fig. 4B) or for fixed mutations (P = 0.08, contingency test; fig. 4C). Thus, we detected no obvious difference in the number or pattern of mutations between rho and rpoB mutations.
We then investigated whether there was evidence of historical contingency at the level of specific mutations. Did rho and rpoB lines tend to accumulate mutations in different sets of genes? We first used a phylogenetic approach: focusing only on mutations that arose during Phase 2 evolution at a frequency of 5% or higher, we calculated a distance matrix and resulting Neighbor-Joining tree from the presence-absence of mutations among populations. We then tested for associations between phylogenetic clustering and the pathways of origin (i.e., rho, rpoB, or control lines evolved from the Phase 1 Ancestor) against the null hypothesis of no associations between pathway and phylogeny. We found significant association between the mutations that arose during Phase 2 and the adaptive history at the level of pathway (Analysis of Similarities: ANOSIM R = 0.139, P = 4 × 10 −4 ; fig. 5A). We also tested the association between mutational patterns and variation across the Phase 2 Founder genotypes, instead of pathways. Here again, the test was significant (ANOSIM R = 0.2398, P = 1 × 10 −4 ). These results are consistent with the idea that the identity of mutations differed among Phase 2 populations based in part on their Founder genotype.
Our second approach relied on Dice's similarity coefficient (DSC) (Dice 1945). Following Card et al. (2021), we calculated DSC at the genic level between all pairs of evolved populations. We based DSC on two sets of mutations: all Phase 2 mutations found at a frequency of 5% or higher (supplementary fig. S4, Supplementary Material online) and Phase 2 fixed mutations found at a frequency of 85% or higher ( fig. 5B). For the former, the average DSC between populations was 0.38, indicating that the evolved populations shared 38% of their mutated genes on average. The mean DSC within a Phase 2 adaptive pathway was 0.39, which was similar to, but significantly different from, the mean DSC between pathways (DSC = 0.37; P = 4.87 × 10 −14 , Wilcoxon test). The same analyses based on fixed mutations had an average DSC 0.10 across all comparisons, an average DSC of 0.06 between populations from different pathways, and a mean DSC of 0.18 within the same pathway ( fig. 5C). The average DSC within and between adaptive pathways again differed significantly (P < 2.2 × 10 −16 , Wilcoxon test), suggesting that fixed (and presumably adaptive) mutations differed among populations in part due to their Phase 2 Founder. We also detected a significant, moderate effect of the adaptive history on DSC across comparison types (η² = 0.078, P < 2.2 × 10-16, Kruskal-Wallis test). Overall, we found that evolved populations from the same adaptive pathway had acquired more similar mutations than populations from differing adaptive pathways (P < 0.01, randomization test). These results further indicate that the fixed mutations that arose in the Phase 2 experiment differed due, in part, to genetic history.
Finally, we sought to identify specific genes that differed in their propensity to house mutations in rho versus rpoB pathways. To do so, we counted the number of populations with and without a mutation in each gene or MBE intergenic region and performed these counts separately for rho and rpoB populations. Using a Fisher's exact test (FET), we identified six genes or intergenic regions that were more frequently mutated in one adaptive pathway but not the other (P < 0.05, table 3). These tests were no longer significant (P > 0.05) after False Discovery Rate (FDR) correction for all but one gene region, likely reflecting low statistical power due to sample size and many (163) FET tests.

Discussion
Evolution is an inherently historical process, but the magnitude and effect of history on adaptation remains somewhat enigmatic. To yield insights into the dynamics of historical contingency, we have performed the second phase of a two-phase evolution experiment. Phase 1 was based on a study that evolved 114 initially identical populations of E. coli to the stressful temperature of 42.2 °C (Tenaillon et al. 2012). These populations evolved primarily by one of two distinct pathways involving mutations in either the RNA polymerase beta subunit gene (rpoB) or the transcription termination factor rho. We chose five clones from each of the two pathways and evolved them for 1,000 generations in a second, low temperature (19.0 °C) environment. At the end of the Phase 2 experiment, we compared the phenotypes (relative fitness, w r ) and genotypes of evolved populations to both their immediate ancestors (Phase 2 Founders) and to the ancestor of the entire experiment (Phase 1 Ancestor; fig. 1).
Based on these data, our first finding is that w r varied significantly among evolved populations both among pathways (rho or rpoB) and among Founder genotypes ( fig. 2), explaining 10% and 23% of the w r variance. Thus, the fitness of Phase 2 populations was contingent upon their Founder fitnesses and genotypes. A more nuanced question is whether there was a predictable pattern to MBE the fitness response. It is reasonable to expect, based on diminishing returns epistasis, that w r among populations vary as a function of the fitness of the Phase 2 Founder, specifically that lower-fitness Founders give rise to populations with more dramatic fitness gains. We have found the expected general trend across rpoB Phase 2 populations ( fig. 3C). Surprisingly, however, this relationship did not hold for the rho lines ( fig. 3C). These contrasting results not only suggest differences among pathways (figs. 3A and B) but raise important questions about what might drive these differences. One must first consider the caveats and limitations of our experimental design. For example, practical considerations limited the number of rho and rpoB Founders; perhaps more or different rho samples would have yielded different results. Moreover, one of the rpoB Founders (rpoB I966S) had a much lower w r than the rest of the Phase 2 Founders, experienced the biggest shift in fitness   . 3C). Unfortunately, we did not have a Phase 2 rho founder for comparison, because the only potential Phase 2 rho founders with similarly low fitness did not survive the 9-day extinction test. Finally, we recognize that the timescale of 1,000 generations likely does not provide a complete view of the fitness landscape; perhaps populations will converge on the same fitness optima with further evolution. The structure of epistasis may drive some of the differences seen among pathways. The ongoing debate about diminishing returns epistasis does not center on whether it exists, because diminishing returns have been found both by testing the effects of individual mutations (Moore et al. 2000;Kryazhimskiy et al. 2009;Perfeito et al. 2014) and by inferring patterns of evolutionary change for experimental evolution data (Chou et al. 2011;Khan et al. 2011;Wang et al. 2018;Bakerlee et al. 2022). Instead, it centers on whether the dynamics of diminishing returns can be predicted based on starting fitness alone or whether diminishing returns is a product of more idiosyncratic epistatic interactions that depend, in part, on the genetic background. As an example of the latter, Card et al. (2019) measured the evolution of antibiotic resistance and found that less-fit genotypes do not always evolve bigger changes in fitness resistance. Like Card et al. (2019), our experiment has not been designed to measure diminishing returns epistasis directly. Nonetheless, our results suggest that fitness evolution is more idiosyncratic than predicted by the global epistasis model, given differences in fitness responses between pathways and genotypes (figs. 2 and 3A-C). Several recent studies have similarly concluded that epistasis is often idiosyncratic (Wei and Zhang 2019;Lyons et al. 2020;Bakerlee et al. 2022).
The next pertinent question is: What might drive these idiosyncrasies? We do not have a complete answer to this question, but we can offer some insights. Previous work has shown that the complete set of 114 high-temperature adapted lines differed substantially in their fitness tradeoffs between 42.2 °C and 19.0 °C (Rodríguez-Verdugo et al. 2014, 2016. A few lines that evolved at 42.2 °C were more fit than their ancestor at 19.0 °C, whereas others were much less fit. Our work further illustrates that Phase 2 populations vary in their trade-off dynamics ( fig. 3D and  E). Furthermore, studies of single mutants have shown that some of the rpoB mutations in this study confer fitness advantages at 42.2 °C, but in contrast, two rho mutations (rho A43T and rho T231A) likely require positively epistatic interactions to become adaptive (Tenaillon et al. 2012;González-González et al. 2017). We suspect that all of these patterns feed into idiosyncratic evolutionary responses. Given, for example, that some genotypes appear to be thermal specialists and others are generalists, one can envision that founding populations with distinct generalist versus specialist mutations have substantially different numbers, directions, and types of potential epistatic interactions across the genome.
Another variable to consider among the Phase 2 Founders is their full genotypic background. As mentioned previously, Tenaillon et al. (2012) illustrated that rho and rpoB lines were positively and negatively associated with their own sets of particular mutations in other genes along the genome. It is possible that other mutations in the genome besides those in rho or rpoB drive the Phase 2 Founder effect. We believe this explanation is unlikely for the rpoB pathway because previous work based on single mutants has shown that most of the rpoB mutations in our Phase 2 Founders cause large shifts in fitness and are responsible for most changes in gene expression observed during hightemperature adaptation (Rodríguez-Verdugo et al. 2016). The case is more nuanced for rho, however, because single rho mutations do not always drive large fitness effects or downstream shifts in gene expression (González-González et al. 2017), suggesting that additional mutations in rho background influence the contingent responses seen here.
A unique feature of our work is that the set of Phase 2 Founders represent pathways that were defined by interactions among distinct set of mutations and genes. Based on the concept of causal dependance (Beatty and Carrera 2011), we predicted that these pathways affect the evolutionary response by shaping the type and identity of future mutations. Although the two pathways do not vary in their number or types of mutations ( fig. 4), there is ample evidence to support our prediction. For example, Phase 2 mutations cluster nonrandomly on a phylogeny ( fig. 5A), suggesting that the set of successful mutations is not independent of the Phase 2 Founding genotype. Similarly, the complement of Phase 2 mutations is more similar within a pathway than between pathways ( fig.  3B). Finally, specific genic and intergenic regions vary in their enrichment for mutations depending on the genetic pathway of their Founders (table 3). [We included intergenic regions because they have been previously implicated as drivers for bacterial adaptation (Khademi et al. 2019).] These patterns hold to some extent for the entire complement of >1,000 mutations, but they are especially clear for the set of 119 fixed mutations ( fig. 5 and table 3). Since fixed mutations are more likely to be adaptive, and also because the number of mutations is unlikely to be limiting in this system, our results show that the identity of adaptive mutations depends on genetic background, likely due to idiosyncratic interactions.
Although our results are not compatible with the global epistasis model (Kryazhimskiy et al. 2014), they do appear to adhere to a modular model of evolution (Tenaillon et al. 2012;Wei and Zhang 2019). Of course, two-phase evolution experiments have inherent biases, because the second phase is always founded by lines already defined by differences in evolutionary outcomes. In this case, we have introduced an additional bias because we have chosen Founders from two distinct pathways. However, this should not inhibit our ability to distinguish between the global or modular models of epistasis in the second phase of evolution. Interestingly, Wei and Zhang (2019) have MBE documented that an emergent property of their modular model is the appearance of global diminishing returns epistasis. That is, epistatic interactions within specific modules can combine to provide the signal of global diminishing returns epistasis.
Our finding-that is, that mutations detected in Phase 2 mutations are associated with Founder pathway-may provide insights about mechanisms of adaptation. We have identified seven genes and intergenic regions that are enriched for Phase 2 mutations in either rpoB or rho populations (table 3). Although only one of the genic regions remained significant after FDR correction, these seven regions constitute a likely set of candidates to drive some our observed genome-wide patterns (based on phylogenetic and similarity analyses) and may therefore yield mechanistic insights. Of the seven, five are more likely to accrue mutations within rpoB populations. In this context, it is worth recalling that rpoB is a component of RNA polymerase, which drives gene expression; as a result, changes in the RNA polymerase have the potential for numerous pleiotropic effects and potential epistatic interactions. At least three of the five enriched regions in rpoB lineages are related to transcriptional function: rpoC, rho, and hepA (table 3). rpoC codes for the beta-prime subunit of RNA polymerase (Trinh et al. 2006;Conrad et al. 2010); the RHO protein terminates RNA polymerase activity; and hepA (which also known as rapA) encodes a transcription factor with ATPase activity that it is an RNA polymerase associated protein (Sukhodolets et al. 2001). Previous work has shown that modifying RNA polymerase is a key feature of adaptation to thermal stress but also that this is a blunt instrument that may cause more phenotypic changes (as measured by gene expression; Rodriguez-Verdugo et al. 2016) than may be necessary to achieve fitness gains. If true, it is reasonable to speculate that adaptation to 19.0 °C from 42.2 °C includes further tuning of RNA polymerase, as reflected by an enrichment of genes related to transcriptional activity like rpoC, rho, and hepA.
Three further features about the enriched regions stand out. First, two enriched regions within rpoB populations (valY and lysV; table 3) encode tRNA synthetases (Andersen et al. 1997;Ruan et al. 2011;Agrawal et al. 2014), suggesting that one additional or alternative route to adaptation is through modifications of translational speed or dynamics. Second, the Phase 2 evolved populations that descended from rho backgrounds were enriched in one gene (ECB_01992) and one intergenic region (ybcW/ ECB_01526) (table 3). Both of these regions have unknown functions, and thus they yield no clues into the molecular mechanisms of 19.0 °C adaptation for rho populations. Finally, it is interesting to speculate about the fact that more enriched regions (5 vs. 2) were found in rpoB vs. rho lines. Previous work has shown that engineered mutations of rho A43T and rho T231A lead to fewer modifications of gene expression than do rpoB mutations I572L, I572N, and I966S (González-González et al. 2017). These observations suggest that some of the rho mutations are "lessconnected" than the rpoB mutations, which could again lead to substantially different dynamics of fitness and epistasis between the two pathways.
To sum, our Phase 2 experiment has shown that the fitness response of evolved populations varies by the pathway and phenotype of their Founders, but Founder fitness is not strongly predictive of fitness gains. More importantly, our work demonstrates that the suite of fixed and presumably adaptive mutations in Phase 2 differs according to their Phase 1 history. Consistent with the concept of causal dependance, our observations illustrate that history has shaped, defined, and perhaps even canalized the adaptive response of Phase 2 populations. Overall, these observations add to a growing literature suggesting that evolution is often contingent on genetic history.

Two-Phase Evolution Experiment Isolate Criteria and Selection
To study evolutionary contingency, we chose ten clones from Tenaillon et al. (2012) (table 1) to evolve at 19.0 °C, which is toward the lower limit of the temperature niche for the REL1206 ancestor (Rodríguez-Verdugo et al. 2014). It is worth noting that the REL1206 ancestor had been evolved for 2,000 generations at 37.0 °C in minimal media prior to the Phase 1 experiment and was therefore preadapted to media and laboratory conditions (Lenski et al. 1991).
Before clones were subjected to the Phase 2 evolution experiment, they were first assessed for survivability. To test survivability, isolates from frozen stock were placed into lysogeny broth (LB) and incubated at 37.0 °C for 1 day to acclimate from frozen conditions (Bennett and Lenski 1993;Lenski and Travisano 1994;Rodríguez-Verdugo et al. 2014). The overnight culture was diluted 1,000-fold in saline, and this dilution was transferred into fresh Davis Minimal (DM) Media supplemented with 25 mg/l of glucose and grown for 1 day at 37.0 °C. Following incubation, 100 μl of the culture was transferred into 9.9 ml of fresh DM media and incubated at 19.0 °C and serially propagated for at least 9 days. Each day, we measured the cell density to determine if extinctions had occurred, by diluting 50 μl of overnight culture into 9.9 ml of Isoton II Diluent (Beckman Coulter) and measuring cell density in volumetric mode on a Multisizer 3 Coulter Counter (Beckman Coulter). An isolate survived if its cell density measurements were maintained over the course of the test whereas allowing for fluctuations of ±1 × 10 6 cells.
Evolution Experiment at 19.0 °C To prepare the isolates for the Phase 2 experiment, the Phase 2 Founders and the Phase 1 Ancestor (REL1206) were grown from frozen stock in 10 ml of LB at 37.0 °C with 120 revolutions per minute (RPM). After 24 h of incubation, the overnight cultures were diluted 10,000-fold and plated onto TA plates and incubated at 37.0 °C. On MBE the next day, single colonies were picked from the plates and inoculated into 10 ml of fresh LB and incubated at 37.0 °C with 120 RPM. The next day, we transferred 100 µl of the bacterial culture into 9.9 ml of fresh DM25 media, which was incubated at 37.0 °C at 120 RPM for 24 h to acclimate to experimental conditions, following common practice (Bennett and Lenski 1993;Lenski and Travisano 1994;Rodríguez-Verdugo et al. 2014). After incubation, we began the Phase 2 evolution experiment by transferring 100 µl of culture into 9.9 ml of fresh DM25 and incubated the tubes at 19.0 °C with 120 RPM and incubated for 24 h.
Each day, the cultures were transferred daily into fresh media via a 100-fold dilution. At regular intervals (at generation 100 and roughly every 200 generations after that), we mixed 800 μl of each line with 800 μl of 80% glycerol to prepare the whole population frozen stocks. We began the experiment in January 2020 but had to pause it after approximately 297 bacterial generations due to the coronavirus disease 2019 pandemic. To restart the experiment, we revived the bacterial populations by transferring 100 μl of thawed glycerol stock into 9.9 ml of fresh DM25 media and continued the experiment until the bacteria had grown for a total of 1,000 generations or 152 days. We note that this restart caused bacterial cultures to experience two carbon sources for 1 day: glycerol and glucose. The use of glycerol may have altered generation time on that day, but we do not expect it had a lasting effect, for two reasons. First, the frozen bacterial stocks were prepared at a final concentration of 40% glycerol that was then diluted 100-fold into fresh DM25 liquid media, such that glycerol constituted a very small proportion of the volume for that day. Second, although the presence of glycerol could affect metabolic function, the selection pressure (temperature) was unaltered.

Measuring Relative Fitness
We performed competition experiments to measure the relative fitness of the Phase 2 evolved lines. We competed the Phase 2 evolved lines against the Phase 1 Ancestor at 42.2 °C and their respective Phase 2 Founders at both 19.0 °C and 42.2 °C. To perform the competitions, we mixed the cells in a single glass culture tube and plated the mixture to count the colonies before and after 24 h of competition. We used the neutral Ara + marker to differentiate between the two lines when plating on tetrazolium-arabinose (TA) plates. To generate Ara + mutants from the Phase 2 Founders for competitions, we followed previously published methods (Lenski et al. 1991).
To validate neutrality, we competed the Ara + mutants against the original Ara-stock using the methods described below. Control competition experiments were performed by competing the Phase 1 Ancestor, E. coli strain B REL1206, against its Ara + mutant, REL1207 (Wiser and Lenski 2015).
To perform competition assays, bacteria from frozen glycerol stocks were revived with a loop into 10 ml of LB and incubated at 37 °C with 120 RPM for 24 h. After incubation, the overnight cultures were vortexed and 100 µl of each were diluted in 9.9 ml of 0.0875% saline solution. From each dilution tube, 100 µl was transferred to 9.9 ml DM25 to incubate at 37.0 °C with 120 RPM for 24 h. Following incubation and in order for the bacteria to acclimate to the experimental temperature, we transferred 100 µl of the overnight cultures into 9.9 ml of DM25 and incubated the tubes at the experimental temperature (19.0 °C or 42.2 °C) with 120 RPM for 24 h (Bennett and Lenski 1993). The next day, we mixed the Ara − and Ara + competitor strains into sterile DM25 media. For competitions at 19.0 °C, we mixed the bacteria 1:1. For competitions at 42.2 °C, we mixed the bacteria 1:1 or we adjusted the ratio to 1:3 if the original ratio resulted in too few colonies (<20) on the plate for either competitor. The mixture was incubated at the experimental temperature of 120 rpm for 24 h. After allowing the cells to compete, we quantified the cell density of each competitor by plating the overnight culture onto TA plates and counting the number of colonies. All competitions were performed in at least triplicate, resulting in roughly 600 competitions.
Using the methods described in Lenski et al. (1991) and Tenaillon et al. (2012), we calculated the relative fitness, w r . The fitness of a Phase 2 evolved line relative to its competitor was estimated by: where E refers to the evolved line and A refers to the ancestral clone, where N E i and N A i represent the initial cell densities of the two competitors, and N E f and N A f represent the final cell densities after 1 day of competition.

DNA Library Preparation and DNA Sequencing
To sequence the evolved populations, we revived populations from ∼10 μl of frozen glycerol stock in 10 ml of DM media supplemented with 1000 mg/l of glucose. The culture tubes were incubated at 19.0 °C with 120 RPM. We extracted from 65 bacterial populations using the Promega Wizard Genomic DNA Purification kit. DNA concentrations were measured with Qubit dsDNA HS Assay kits. We prepared our DNA sequencing libraries with the Illumina Nextera DNA Flex Library Preparation kit. The libraries were multiplexed and sequenced using the Illumina NovaSeq on an S4 flow cell to generate 100 bp paired-end reads at UC Irvine's Genomics High-Throughput Facility (https://ghtf.biochem.uci.edu). Sequencing read quality was assessed with FastQC v. 0.11.9 (http://www. bioinformatics.babraham.ac.uk/projects/fastqc), trimmed with fastp v. 0.23.2 (Chen et al. 2018), and visualized with MultiQC v. 1.9 (Ewels et al. 2016). Each population had 25,000,000 sequencing reads on average (min = 8,600,000 reads, max = 35,000,000 reads), resulting in a minimum of >150× coverage per population.

Variant Detection
We detected mutations and their respective frequencies in each evolved Phase 2 population using breseq v. 0.35.5 (Deatherage and Barrick 2014). We performed the breseq analysis in polymorphism mode with two different reference genomes. First, we performed breseq analysis using E. coli strain B REL606 as the reference genome. This E. coli strain differs from the Phase 1 Ancestor, REL1206, in seven positions (topA, spoT K662I, glmU/atpC, pykF, yeiB, fimA, and the rbs operon) that were excluded from our analysis (Barrick et al. 2009;Tenaillon et al. 2012). We performed a first round of variant detection using breseq in polymorphism mode on all of evolved populations relative to E. coli strain B REL606.
Following this first step of the analysis, we generated a new, mutated reference sequence to represent each Phase 2 Founder using the gdtools APPLY command in breseq using the sequencing data available in Tenaillon et al. (2012). We then ran the breseq analysis again with respect to the Phase 2 Founder using the mutated references to verify mutation predictions, as described in Deatherage and Barrick (2014). Using gdtools available through breseq, we compiled the mutation information into readable tables and as an alignment file in PHYLIP format. A phylogeny was constructed using IQ-tree and the PHYLIP alignment as input (Nguyen et al. 2015).

Statistical Analyses
All statistical analyses were performed in R v 4.0.2 (R Core Team 2019). For relative fitness results and statistical analysis, we first assessed the normality of the data using the Shapiro-Wilk test and the variance with Levene's test available through R. To perform ANOVA, ANCOVA, and Kruskal-Wallis analyses, the R package rstatix v 0.7.0 was used (Kassambara 2023). Model II regression analyses were also carried out in R using the package lmodel2 v 1.7-3 (Legendre 2018). To statistically test for associations between the mutation patterns observed in Phase 2 and their initial adaptive pathway, we first built a distance matrix from the presence and absence matrix of Phase 2 mutations in R. Using the vegan package v 2.5-7 in R, we directly tested for associations between the distance matrix of Phase 2 mutations and the adaptive pathway or mutated codon with ANOSIM (Oksanen et al. 2020). We also built a Neighbor-Joining (NJ) tree based on the presenceabsence matrix of accessory genes. To do so, we first calculated the Euclidean distances from the presence-absence matrix of the accessory genes using the dist function in R. We then built the NJ tree from the Euclidean distances using the ape package in R (Paradis and Schliep 2019). We calculated DSC for each pair of Phase 2 evolved populations using the vegan package. Dice's coefficient of similarity was calculated using the same distance matrix used to build the NJ tree described above, as well as on a distance matrix containing only the fixed mutations. In order to statistically analyze Dice's coefficient of similarity across our groups, we employed a randomization test following previously established methods (Card et al. 2021;Debray et al. 2022). First, we calculated the average of the coefficients for all comparisons within an adaptive pathway and for all comparisons between pathways. The difference between these values was calculated and compared against simulated means. We generated simulated means to serve as a null hypothesis by randomly shuffling the similarity data with 5,000 iterations. We then calculated significance by measuring the proportion of permutations in the expected distribution with a specificity statistic value greater than or equal to the observed value (Card et al. 2021). To identify genes or intergenic regions that were more frequently mutated in populations descended from one adaptive history but not the other, we built a 2 × 2 FET for each mutated gene or intergenic region in R, for a total of 163 contingency tests.

Supplementary Material
Supplementary data are available at Molecular Biology and Evolution online.