Strong selective environments determine evolutionary outcome in time‐dependent fitness seascapes

Abstract The impact of fitness landscape features on evolutionary outcomes has attracted considerable interest in recent decades. However, evolution often occurs under time‐dependent selection in so‐called fitness seascapes where the landscape is under flux. Fitness seascapes are an inherent feature of natural environments, where the landscape changes owing both to the intrinsic fitness consequences of previous adaptations and extrinsic changes in selected traits caused by new environments. The complexity of such seascapes may curb the predictability of evolution. However, empirical efforts to test this question using a comprehensive set of regimes are lacking. Here, we employed an in vitro microbial model system to investigate differences in evolutionary outcomes between time‐invariant and time‐dependent environments, including all possible temporal permutations, with three subinhibitory antimicrobials and a viral parasite (phage) as selective agents. Expectedly, time‐invariant environments caused stronger directional selection for resistances compared to time‐dependent environments. Intriguingly, however, multidrug resistance outcomes in both cases were largely driven by two strong selective agents (rifampicin and phage) out of four agents in total. These agents either caused cross‐resistance or obscured the phenotypic effect of other resistance mutations, modulating the evolutionary outcome overall in time‐invariant environments and as a function of exposure epoch in time‐dependent environments. This suggests that identifying strong selective agents and their pleiotropic effects is critical for predicting evolution in fitness seascapes, with ramifications for evolutionarily informed strategies to mitigate drug resistance evolution.

and 200 µM of dNTPs in a final volume of 50 µL of 1 × Standard Taq Reaction Buffer. The cycling conditions were as follows: 95 ºC for 5 min to enhance cell lysis, and 30 cycles of 95 ºC for 30 s, 56 ºC for 30 s and 72 ºC for 1 min, with a final extension at 72 ºC for 5 min. The PCR products were Sanger sequenced by a third party (Institute of Biotechnology, University of Helsinki, Finland) using in-house protocols. The sequences were assembled from Sanger sequencing chromatograms using Pregap4 and Gap4 in the Staden Package (Staden et al. 2000), and aligned with ClustalW using MEGA7 (Kumar et al. 2016) (Fig. S1).

Competition assay for ancestral strains: E. coli REL606 and REL607-like Ara + reversion mutant
To ascertain that there was no initial difference in fitness between the two ancestral strains, a competition assay was carried out following a previously established protocol (barricklab.org/twiki/bin/view/Lab/ProceduresLongTermCompetitions, accessed 2018-05-16). For the assay, the strains were revived by culturing 2 µL of freeze-stored clonal culture in 5 mL LB liquid medium overnight at 37 ºC with constant rotation at 120 r.p.m. Subsequent culturing steps were performed in DM1000 medium. Before starting the competition assay, the strains were mixed together at 200-fold dilutions (100-fold overall dilution in cell number), and this mix was immediately dilution plated on TA agar to determine the initial frequencies of the strains. Subsequently, 50 µL of each of the strains was added to 10 mL medium, with 12 replicates, followed by culturing at 37 ºC, and serially diluting (1 % volume) each day for a total of three days. This allows relative fitness to be determined with a precision of ± 1 % (95 % confidence intervals). Upon completion of the assay, the cultures were immediately dilution plated on TA agar to determine the final frequencies of the strains. agar, DF the dilution factor of all transfers combined, and i and f the initial and final time points, respectively. A difference between the logarithm of the Malthusian parameters (ratio ≠ 1) for the two strains (paired for each replicate) indicates a difference in fitness. The result of the competition assay (one of 12 replicates was lost due to technical error, i.e. N = 11) supports the absence of a difference in fitness between the two ancestral strains (Fig. S2). 5000 bootstrap samples were taken; the confidence interval is bias-corrected and accelerated. The Pvalue reported is the likelihood of observing the effect size if the null hypothesis of zero difference is true.

Details on serial passage experiment
The experiment was started by culturing the REL606 and REL607-like strains in DM1000 for 24 h to obtain 2 × 10 9 cells mL −1 . One colony from each plated strain was added to 2 × 5 mL of DM1000, and from this solution, approximately 10 6 bacterial cells (50 µL) of the corresponding strain were pipetted to the wells of the deep well plates containing 500 µL of DM1000 with 0.5 × MIC of the appropriate antimicrobial. In addition, 5 × 10 5 plaque forming units (PFU) of the T4 phage (13 µL from phage stock containing 3.9 × 10 9 PFU mL −1 ) were pipetted to phage treatment wells. Deep well plates were cultured at 37 ºC with constant shaking at approx. 120 r.p.m. Serial passage was performed every 24 h by transferring 10 µL (2 % v/v) from each well to a new well containing fresh medium using a multichannel pipette. To store populations in suspended animation, after the first 48 h and subsequently every 96 h in the experiment, all deep well plates were freeze-stored at -20 ºC by mixing in 250 µL of sterile 85 % glycerol in each well and covering with a foil seal.

Determining minimum inhibitory concentrations for antimicrobials for ancestral E. coli strain
The MIC values for the experimental antimicrobials were determined for the REL606 strain in the experimental conditions using the microdilution method. Instead of 2-fold dilutions typically used, the experiment was performed at even concentration intervals in a total interval determined based on expected MIC values: for nalidixic acid, 1-5 µg mL −1 (expected MIC of µg mL −1 for REL606 based on (Harmand et al. 2018); for rifampicin, 1-5 µg mL −1 (expected MIC of 2.5 µg mL −1 for E. coli B (Rodriguez-Verdugo et al. 2013)); for spectinomycin, 10-20 µg mL −1 (expected value of 10-20 µg mL −1 for E. coli B (Anderson 1969)). The experiment was performed by adding 10 6 cells (1/2000 dilution of overnight culture of REL606 in DM1000 medium) to honeycomb wells containing DM1000 medium, with each well containing a different concentration of the antimicrobial in question according to the gradients specified above. This cell number is expected to be lower than the spontaneous frequency of resistance mutations, indicating that any observed growth represents the intrinsic resistance level of the ancestral strain rather than that of derived mutants emerging during the test. Culturing was performed for 24 h at 37 ºC using the Bioscreen C well-plate reader (Labsystems, Helsinki, Finland) that measures optical density (OD) at 420-580 nm with a wideband filter at 5 min intervals.
The MIC value was interpreted as the smallest concentration preventing growth, which was 1.0, 2.6 and 8.0 µg mL −1 (experimental concentrations of 0.5 × MIC thereby being 0.5, 1.3 and 4.0 µg mL −1 ) for nalidixic acid, rifampicin, and spectinomycin, respectively (Fig. S3). The dose response curves in Fig. S3 show that the antimicrobial level where the carrying capacity begins to decrease is lower than the concentration (0.5 × MIC) used in the experiment. This indicates that differences in effective population size caused by the antimicrobials do not account for differences in evolutionary dynamics between the treatments in the absence of evolution altering the population size response to the antimicrobials.

Measuring antimicrobial resistance phenotypes
To quantify the evolution of antimicrobial resistance over time, every 2 days in the 48-day evolutionary experiment, 2 % (10 µL) of the old culture was transferred, using a 96-pin replicator, to a large petri dish containing lysogeny broth (LB) agar supplemented with one of the three antimicrobials at a selective concentration (> 1 × MIC), including an antimicrobial-free plate to control for bacterial extinction. After culturing for 24 h at 37 ºC, growth on the plates was quantified on a binary scale (0 = no growth; 1 = growth). Therefore, for this metric, growth indicates the detectability of antimicrobial resistance in a minimum of 2 % of the bacterial population. Notably, we cross-phenotyped all populations from all treatments and time points against all antimicrobials. Since the temporal tracking of resistance evolution employed a coarse measure diverging from experimental conditions (solid medium vs. liquid medium; higher antimicrobial level compared to selective concentration of 0.5 × MIC; phenotype based on fraction of heterogeneous population), we isolated individual clonal lineages from the experimental end-point and phenotyped them more precisely. Clones were isolated by streak-plating inoculum from the population on an LB agar plate, culturing overnight, and picking an individual colony, followed by propagation in liquid medium. The clones were phenotyped by culturing a small initial number of cells for 24 h in DM1000 liquid medium without antimicrobials or supplemented with one of the three antimicrobials at a range of multiplicities of the MIC value of the ancestral E. coli strain: 1, 2, 10, or 100 × MIC. This allowed obtaining quantitative growth and antimicrobial resistance phenotypes (in units of optical density at 600 nm, OD600) in the experimental conditions for isogenic lineages likely to represent the dominant E. coli genotype in each of the 900/928 populations that escaped extinction during the 48-day serial passage experiment.
Before downstream analyses, the OD data was curated by removing background OD values specific to 96-well plate location, drug type, and drug concentration (Fig. S4), and following this, by removing (transforming to zero growth) likely false positive values based on belonging to a distinct distribution close to zero (Fig. S5). The cut-off used to remove false positives was 3 times the standard deviation of the noise distribution close to zero, representing a demarcation region between the noise distribution and true positive value distribution. In later analyses, resistance was defined for the clones as the capacity to grow at antimicrobial concentrations exceeding the MIC of the ancestral strain, as the experimental concentration of 0.5 × MIC selected mostly for low-level resistance (Fig. S6).

Measuring phage presence and resistance phenotypes
We also measured the loss of phage from phage-containing populations during the experiment using the population level pin replicator system on agar plates containing a soft agar lawn of the ancestral E. coli strain, where the formation of a plaque indicates presence of phage. Moreover, we quantified the phage resistance phenotype of the end-point clones using a similar setup, with a lawn of the ancestral T4 phage instead of the bacterium. Bacterial growth on the phage lawn was taken to indicate phage resistance, which was therefore characterized as a binary trait (susceptible/resistant). Detailed phenotyping protocols are documented below.

Plaque assay and generation of T4 phage stock
The T4 phage stock for the experiment was generated based on previously established protocols (Lenski lab protocol: lenski.mmg.msu.edu/ecoli/phagepharm.html, accessed 2018-05-15; Barrick lab protocol: barricklab.org/twiki/bin/view/Lab/ProceduresPhageFarming, accessed 2018-05-18). First, a dilution series of an existing T4 stock was prepared in LB liquid medium, and 100 µL was combined with 100 µL of exponentially growing host strain (REL606). Following incubation at 37 ºC for 30 min, the mixture was added to 4 mL of soft LB agar (7 g agar per litre) kept at 55-60 ºC, mixed and poured on top of an LB agar plate. The plates were cultured for 24 h at 37 ºC, generating phage plaques at dilutions where the plaques are not touching. Subsequently, 1 mL of overnight host culture was transferred to 4 mL of LB medium in a 50 mL flask, and a toothpick was stuck into an isolated plaque and dropped into the host culture, which was kept at 37 ºC with constant shaking at 150 r.p.m. After 6-8 h causing complete lysis of the host, 10 drops of chloroform were added and mixed to kill the remaining bacteria, and cell debris was removed by spinning down and keeping the supernatant. The phage was then titred as described above to determine phage density in the stock (as plaque forming units, or PFU, i.e. infectious phage particles). The lysate was stored with glycerol at -80 ºC.

Details on phenotyping experimental populations over time
To track the development of antimicrobial resistance over time, the following selective concentrations were chosen for each antimicrobial for use in LB agar plates based on literature: 2 and 20 µg mL −1 for nalidixic acid, 50 µg mL −1 for rifampicin, and 50 and 75 µg mL −1 for spectinomycin. After each 24 h growth cycle in the experiment, the populations were transferred to selective plates, as well as plates without antimicrobials to control for the presence of bacteria, by pin-replicating the whole deep well plate onto a large petri dish. To assess the presence of phage, 200 µL of overnight culture of the ancestral REL606 strain into was added to a sterile 15 mL Eppendorf tube, followed by adding 15 mL of soft agar tempered to 55 ºC in a water bath, vortexing, pouring on large petri dish containing LB agar, and pin-replication of the 24 h culture from the experiment upon solidification of soft agar. Plates were cultured for 24 h at 37 ºC. Each plate was documented by photography, and the results were interpreted as negative (0), weak (1), or strong (2) bacterial growth (antimicrobial resistance) or plaque formation (phage presence), subsequently converted to binary data (2 converted to 1).

Details on phenotyping clones isolated from experimental end-point
Clones were isolated from the experimental end-point by streaking an inoculum from surviving populations (N = 900 out of 928 populations in total) on LB plates. Following culture (24 h / 37 ºC), a single colony was selected and transferred to a 96 well plate containing 200 µL of DM1000. Following culture (24 h / 37 ºC / 120 r.p.m.), the clonal lines were freeze-stored with 85 % glycerol at -80 ºC. To perform subsequent tests, the master well plates containing the original clones were replicated on 96 well plates containing 200 µL of DM1000. The test well plates were cultured and freeze-stored as above. To obtain quantitative growth traits (optical density at 600 nm) for different antimicrobial concentrations, a test well plate was replicated on another 96 well plate containing 200 µL of DM1000, and cultured for 24 h at 37 ºC / 120 r.p.m. This results in a cell density of 2 × 10 9 cells mL -1 . The clonal cultures were diluted into deep well plates to obtain a total dilution of 1:2000, corresponding to 1 × 10 6 cells mL −1 assuming a density of 2 × 10 9 cells mL −1 . This cell count is expected to be lower than the frequency of spontaneous antimicrobial resistance mutations, indicating that any observed growth represents the (potentially evolved) resistance level of the clonal line rather than that of derived mutants emerging during the test. Approx. 2000 cells (2 µL) of each diluted clone was pipetted on another 96 well plate containing 300 µL of DM1000 with each of the experimental antimicrobials at the following concentrations: 0, 1, 2, 10, or 100 × MIC (respectively: 0, 1.0, 2.0, 10, and 100 µg mL −1 for nalidixic acid; 0, 2.6, 5.2, 26, and 260 µg mL −1 for rifampicin; and 0, 8.0, 14, 80, and 800 µg mL −1 for spectinomycin. Following culture for 24 h at 37 ºC / 120 r.p.m., the cell density reached by each clone was quantified by measuring OD at 600 nm wavelength (Tecan Infinite M200) with a bandwidth of 9 nm and 25 measurements per well. The REL606 and REL607-like strains and blank wells containing the medium alone supplemented with the different antimicrobials (only rifampicin was shown to alter the background OD level) were used as controls in each well plate. To determine the phage resistance phenotype of the clones from the phage treatments, a work plate containing REL606 and REL607-like strains as controls was pin-replicated on an LB plate containing a top layer of soft agar with the ancestral T4 phage (see detailed protocol for plaque assay above). Following culture for 24 h at 37 ºC, phage resistance was quantified for each clone as a binary trait based on the absence (0) or presence (1) of bacterial growth on the phage lawn.

Details on extraction of DNA from clones isolated from experimental end-point
DNA extraction was performed from 1 mL of overnight culture with the Qiaqen DNeasy 96 Blood & Tissue kit using a custom protocol. To release nucleic acids, the cells were transferred on S plates and spun down for 30 min at 6200 r.p.m., followed by addition of 200 µL of ATL and proteinase K solution (prepared by mixing 36 mL of ATL and 4 mL of proteinase K) to each well, vortexing for 15 s with plastic cover attached, and briefly spinning down at 3000 r.p.m. The suspension was incubated for 30 min at 56 ºC, with vortexing every 10 min. Subsequently, 410 µL of AL solution was added to each well, followed by vortexing for 15 s, and spin down at 3000 r.p.m. as described above. The whole sample (approx. 620 µL) was transferred to the extraction membrane columns, followed by spinning down for 10 min at 6000 r.p.m. in room temperature (RT) with plastic cover attached and discarding the flow-through. Subsequently, 500 µL of AW1 was added, followed by spinning down at 6000 r.p.m. in RT and discarding the flow-through. The plastic cover was removed, followed by adding 500 µL of AW2, spinning down for 15 min at 6000 r.p.m. in RT, discarding the flow-through, and centrifuging empty wells for 2 min at 6000 r.p.m. in RT to remove ethanol residues. After the ethanol was allowed to evaporate for 2 min in RT, membrane columns were transferred to elution plates.
Elution was performed in DNA-free grade water by incubating in 40 µL of elute for 2 min in RT, centrifuging for 2 min at 6000 r.p.m., adding another 40 µL of sterile DNA-free grade water and centrifuge immediately for 2 min at 6000 r.p.m. The extracted DNA was stored at -20 ºC.

Genome alignment, variant calling, and annotation
Quality controlled FASTQ files were aligned to the reference genome (NCBI Reference Sequence NC 012967, assembly ASM1798v1) with Bowtie 2 v2.3.4 (Langmead and Salzberg 2012) using default settings. SAMtools v1.4 (Li et al. 2009) was used to convert thus obtained SAM files to BAM files, and to sort and index BAM files. Picard v2.18.10 (http://broadinstitute.github.io/picard) was used to mark duplicates and add read group information to BAM files, and following these steps, to compute genome

Details on machine learning analysis
Before analyses, control treatments were excluded and rare features were removed based on min. 90 % of values being zero, and the data was standardized by converting each value into a Z-score (subtracting each sample's mean and dividing by the sample's standard deviation). The package splitstackshape (Mahto 2019) was used to assign 80 % of each antimicrobial history into a training set and the rest 20 % into a test set. This was followed by random forest classification using the function randomForest implementing the Breiman's random forest algorithm, with the options importance = TRUE and proximities = TRUE. The model was evaluated on the test set using the predict function, and the observed and predicted data were used to produce a confusion matrix with the confusionMatrix function in the caret package (Kuhn 2016). To avoid prediction bias due to a small test set size, this procedure was iterated 100 times, and the resulting confusion matrices were averaged over the iterations to produce one aggregate confusion matrix. The predictability of each estimated factor was interpreted as the weight of the diagonal (observed value matches predicted value) in the confusion matrix.

Details on potential molecular targets of resistance
Nalidixic acid resistance: Nalidixic acid exposure was associated with a small number of nonsynonymous mutations overall (median of 0 mutations in time-invariant single agent environment in the absence of phage, equivalent to control environment) as well as a relatively small number of recurrent mutational targets occurring only in a small proportion of the isolates (Fig. 2G). This is consistent with major clones lacking resistance in time-invariant nalidixic acid environment in the absence of phage (Fig. 2C). These targets include acrR (encoding HTH-type transcriptional regulator), rfaQ (lipopolysaccharide core heptosyltransferase) and ECB RS03400 (putative phosphoglucomutase). Among these genes, acrR has been previously implicated in quinolone resistance (Schneiders et al. 2003), while the product function (LPS biosynthesis) and previous findings suggest that rfaQ (Girgis et al. 2009) and phosphoglucomutase (Paterson et al. 2009;Correia et al. 2014) may be associated with resistance to both quinolone and phage.
Rifampicin resistance: The vast majority of clones from the time-invariant single agent environment had nonsynonymous mutations (median of 1.5 mutations in the absence of phage) in the gene rpoB (β subunit of RNA polymerase) known to produce rifampicin resistance (Mariam et al. 2004), and the presence of mutations in this gene was almost exclusively associated with a resistance phenotype (Fig. 2G). In addition, five other genes (galU, infB, marR, mreC and mrdB) were recurrently mutated only in the presence of rifampicin. Among these, marR is related to drug efflux (Grove 2013); mreC and mrdB are both related to cell shape, with mreC mutations having been previously associated with rifampicin exposure (Gliniewicz et al. 2015); and mutations in the essential gene infB (translation initiation factor 2) have been previously implicated in rifampicin resistance (Huseby et al. 2020) and compensating for the fitness cost of antimicrobial resistance mutations (Zorzet et al. 2010). The gene galU is involved in LPS biosynthesis and has been previously implicated in phage resistance (Le et al. 2014), and consistent with this, was only hit in the presence of phage in this study.
Spectinomycin resistance: The time-invariant spectinomycin (aminoglycoside class) environment, either in the presence or absence of phage, led to higher mean spectinomycin resistance prevalence among end-point clones compared to the antimicrobial-free control and time-invariant combination environments, although resistance increased markedly only in the time-dependent environment in the absence of phage ( Fig. 2E; antimicrobial regime, P < 0.001; for full results, see Table S5). We were unable to obtain robust population-level time series data for spectinomycin, as the majority of populations in all test conditions (time points and spectinomycin concentrations) gave a positive result with the pin replication on agar plate method, masking any potential differences between experimental treatments. There are a number of factors that could enable population growth in our selective conditions, resulting in a positive resistance signal despite all or most of the cells remaining sensitive.
Resistance may occur only in a subset of the population owing to a low selection coefficient, reversible amplifications Nicoloff et al. 2019), or plastic gene regulatory changes induced by antimicrobial stress (e.g. SOS or stringent response) (Poole 2012;Strugeon et al. 2016). Three genes, nadR (transcriptional regulator), trkH (potassium uptake protein) and ECB RS09520 (putative carboxyl-terminal processing protease), were recurrently mutated specifically in the presence of spectinomycin (median of 1 nonsynonymous mutation in time-invariant single agent environment), with only trkH reaching high (close to 0.5) frequencies (Fig. 2G). Among these, trkH has been previously implicated in aminoglycoside resistance (Lazar et al. 2013;Oz et al. 2014), while nadR (Zhang et al. 2018) and proteases (Culp and Wright 2017) can, among other functions, be related to the bacterial stress response. The observations that resistance development was difficult to distinguish in the timeseries data and that two of the three potential molecular resistance mechanisms detected were linked to generic stress response rather than being specific to the antibiotic agent could be indications of weak selection pressure for resistance in the time-invariant spectinomycin environment. As much higher resistance levels occur in a subset of the time-dependent environments, with a shorter spectinomycin selective window, these are expected to arise from the other agents (see results section "Pleiotropic and fitness effects of strong selective agents, rifampicin and phage, modulate evolutionary outcome" in main text).
Phage resistance: Phage exposure caused recurrent mutations in several genomic targets almost exclusively associated with phage resistant phenotypes, with most targets exhibiting low to moderate (max. 0.5) prevalence among isolates (Fig. 2G). This indicates that resistance to the phage T4 has a wide target in E. coli B. This is consistent with earlier studies showing that membrane modifications preventing phage adsorption represent a highly common resistance mechanism of bacteria against virulent phages and can typically be achieved by mutations in a number of genes affecting membrane structure and components (Labrie et al. 2010). Moreover, the phage caused an increase of one nonsynonymous mutation in the median mutation count of the clones, suggesting that individual mutations rather than several mutations in combination were required for phage resistance. The potential phage resistance targets discovered in this study encompass 14 genes: acrR, asmA, fabR, galU, infB, lpcA, mscM, rfaQ, ECB RS0200, ECB RS03400, ECB RS03925, ECB RS05770, ECB RS09520, and ECB RS18465 (Fig.   2G). Six among these (acrR, galU, infB, rfaQ, ECB RS03400 and ECB RS09520) were also associated with resistance to a particular antimicrobial compound and have been discussed above. Among the eight remaining genes specifically associated with phage exposure and phage resistance phenotypes, mutations in four genes (asmA, fabR, mscM and ECB RS18465) were particularly prevalent, with mutations in mscM being prevalent across antimicrobial compound environments but absent from the antimicrobial-free control environment. The product of asmA is involved in the assembly of outer membrane proteins and has been implicated in phage resistance (Misra and Miao 1995). The product of medium with 1000 mg L −1 glucose (DM1000) for culturing strains in liquid medium during the competition assay, serial passage experiment, and clone measurements. Culture conditions: 24 h at 37 º C. The medium is prepared by adding dH2O to final volume and autoclaving, followed by adding 0.5 mL each of the following stock solutions: 10 % (w/v) magnesium sulfate MgSO4 (separately autoclaved stock) and 0.2 % (w/v) thiamine (vitamin B1; filter sterilized), and 5 mL of 10 % glucose (separately autoclaved stock). Sources: Levin et al. (1977)   The water is split into three parts, salts, agar, and sugar, which are autoclaved separately and combined, followed by the addition of 1 mL of each of the following stock solutions: 10 % (w/v)  Figure 2, Figure 3, Figure 4, and Figure 5. Recipe for tetrazolium and arabinose (TA) agar plates for distinguishing between Ara − (red/purple) and Ara + (white/pink) strains. Culture conditions: 24 h at 37 ºC. The medium is prepared by splitting the water, autoclaving the arabinose separately, combining the two parts after autoclaving, adding TTC from a sterile stock solution, and mixing. Sources: Carlton and Brown (1981) Figure 2, Figure 3, Figure 4, and Figure 5. ANOVA table for generalized least squares model on temporal nalidixic acid resistance dynamics at population level.