Temperature-induced prophage dictates evolution of virulence in bacteria

2 Environmental changes can cause strong cascading effects in species communities through altered species interactions. The highly specialized interactions 4 arising from the co-evolution of hosts and parasites such as bacteria and phages as well as their short generation times can rapidly lead to considerable evolu6 tionary responses to abiotic changes, with potential large-scale ramifications to other trophic levels. Here we report how the temperature-dependent induction of 8 a temperate phage in the bacterium Serratia marcescens resulted in both genetic and epigenetic changes in the bacterial host after experimental evolution under 10 different temperatures, and subsequently altered the virulence of the bacterium itself in an insect host. The majority of the genetic and epigenetic changes were 12 associated with the indirect effects of the phage on the bacterium rather than with the experimentally imposed thermal environments. Given the abundance 14 of temperate phages in bacterial genomes, the sensitivity of the onset of their lytic cycle to environmental conditions, and the predominance of environmental 16 changes due to climate change, our results warrant attention as a cautionary example of the dangers of predicting environmental effects on species without 18 considering complex biotic interactions.


22
Understanding the consequences of environmental changes on ecosystems is of high priority as climate change alters environmental conditions (Chevin et al., 2010;Kristensen 24 et al., 2020). Abiotic environmental selection pressures can affect species directly or indirectly via their ecological interactions with other more sensitive species. Under-26 standing those indirect effects is critical as they can amplify through ecosystems (Zarnetske et al., 2012) and requires studying multi-species systems rather than focusing 28 on single species. In this respect, phage-bacteria systems are both relevant for natural ecosystems and practical for laboratory experiments. Phages are ubiquitous pathogens 30 that can drive the abundance and composition of bacterial communities (Suttle, 2005).
Many phages exist as inducible prophages present in bacterial genomes and environ-32 mental perturbations can trigger their induction (Canchaya et al., 2003;Williamson et al., 2002;Brum et al., 2016;Howard-Varona et al., 2017). Since bacteria have a 34 key role in nutrient and energy cycles, indirect consequences of prophage induction can range from modulation of biogeochemical cycles (Suttle, 2005) to life and death of 36 multicellular hosts of pathogenic bacteria (Buck and Ripple, 2017).
We used an evolution experiment and single-molecule real-time (SMRT) sequencing 38 to determine the effects of temperature changes on the evolution of a bacterium (Serratia marcescens) carrying a temperature-sensitive prophage. We describe a novel trophic 40 cascade where prophage induction directed the evolution of its bacterial host, and show that prophage-driven evolution cascaded across trophic levels by having strong effects Figure 1: Effect of evolutionary treatment and assay temperatures on the estimated induction of prophage PP4. Assays lasted two days and assay temperatures are given as day1/day2. A) posteriors of the model-estimated mean for each treatment/assay combination. Points are estimated phage release rates for each of the 29 sequenced clones. B) estimates of the evolutionary treatment effects and C) estimates of the assay effects, with the reference strain used as a reference point. Posteriors are shown as median and 95% credible interval. One-sided Bayesian p-values for pairwise comparisons denoted by * (p < 0.05) and *** (p < 0.001).  several candidate strains, including the reference strain, but no plaque was observed in any of the conditions tested. Taken together, those results confirm that the PP4 76 prophage is induced in S. marcescens cells and released as phage particles at a low rate, as detected by qPCR and confirmed by TEM, and that the induction and/or 78 infectivity rates of this phage are low enough that no clearing is observed in our plaque assays. Interestingly, the immunity repressor protein located immediately upstream of The induction of prophage PP4 as estimated by our method was more pronounced in the reference strain at 24 • C and 31 • C rather than at 38 • C. This environmental 88 sensitivity suggested that cooler environments could act as a strong selective force in S. marcescens as evolution in cooler environments could potentially trigger counter-90 adaptations in bacteria to reduce the fitness loss due to prophage induction and cell lysis (Canchaya et al., 2003), while evolution in warmer environments could on the 92 contrary release selective pressure against such fitness loss.
To test if environmental changes could affect the relationship between prophage 94 PP4 and its bacterial host, we used bacterial strains from an earlier experimental evolution study where temperature had been manipulated (Ketola et al., 2013): cultures of 96 S. marcescens had been left to evolve under either (i) constant moderate temperature at 31 • C, (ii) constant high temperature at 38 • C, or (iii) daily fluctuating temperature 98 between 24 and 38 • C (Supplementary Figure S2). We selected several independently evolved clones from each evolutionary treatment (n = 10, 8, and 10, respectively) 100 and measured prophage PP4 induction with two-day thermal assays ( Figure 1A). The thermal assays were (daily temperatures given for first/second day) 24/24 • C, 24/38 • C, 102 31/31 • C, 38/24 • C and 38/38 • C, thus enabling to test both the effect of mean temperature and of temperature fluctuations on prophage induction. 104 As hypothesized, we found an evolutionary adjustment of the prophage induction in cooler environments (evolutionary treatment effect, Figure 1B). The strains evolved 106 at 31 • C and at 24-38 • C released 19% and 24% less phages, respectively, than strains evolved at 38 • C. Mean patterns of PP4 induction did not differ significantly between 108 strains that had evolved at lower mean temperature (31 • C versus 24-38 • C). The main driver of prophage induction in our assays was the temperature of the final day of the 110 assay, rather than whether a temperature change was experienced during the assay (assay temperatures effect, Figure 1C). Ending an assay at 31 • C induced about three 112 times more phages than ending an assay at 24 • C, and ending an assay at 24 • C induced about ten times more phages than ending an assay at 38 • C. Higher prophage induc-114 tion rates at lower assay temperature are likely to have an adverse effect on bacterial fitness in laboratory culture conditions via cell lysis. The evolution of lower prophage 116 induction rates in the evolutionary treatments where those temperatures occurred is consistent with such a temperature-dependent effect of PP4 on S. marcescens fitness.

Induced evolutionary changes have indirect consequences for bacteria's hosts
After confirming the evolutionary changes in the interaction between prophage PP4 120 and its bacterial host, we explored how these environmentally-triggered changes might cascade to another trophic level: a host of the bacterium itself. We estimated the viru-122 lence of the experimentally evolved strains in an insect host by measuring the survival time of waxmoth larvae (Galleria mellonella) infected by an injection of 5 µl of bacte-124 rial cultures and placed into two assay environments: 24 • C and 31 • C (Supplementary Figure S4). We did not use 38 • C as this temperature would have been lethal to wax-126 moth larvae. Control larvae were injected with sterile water and sterile medium. We estimated the relative virulences of sequenced strains from larvae survival data with a

128
Cox proportional hazards mixed model which controlled for larva body mass and for the initial density of the bacterial cultures and which allowed for different variances 130 across evolutionary treatments ( Figure 3A). Those estimates revealed that the average virulence of clones evolved at high temperature (38 • C) tended to be higher than for 132 clones evolved at lower mean temperature when larvae were incubated at 24 • C (38 • C versus 24-38 • C, one-sided Bayesian p-value = 0.057; 38 • C versus 31 • C, p = 0.066; 134 Figure 3B). However, when larvae were incubated at 31 • C, these differences disappeared while clones evolved at 31 • C were more virulent than those evolved at 24-38 • C 136 (p < 0.01). To confirm those tentative results obtained from the 29 sequenced clones, we utilized a much larger pool of evolved clones which confirmed that overall clones 138 evolved at 38 • C had indeed a higher virulence than the others when assayed at room temperature (p < 0.01 for comparisons of 38 • C clones with both 24-38 • C and 31 • C 140 Figure 3: Effect of evolutionary treatment on strains virulence in waxmoth larvae at two incubation temperatures. A) relative virulence of individual clones, measured as relative hazards estimated from a Bayesian implementation of a Cox proportionalhazards model. All virulence estimates are relative to the virulence of the reference strain in incubation at 24 • C (denoted by a broken horizontal line) and are corrected for the effects of injection batch, larval body mass and optical density of injected cultures. B) Mean relative virulence per evolutionary treatment and per incubation temperature as estimated by the model (exp(µ evo )) (n = 29 sequenced clones). C) Confirmatory results from a similar virulence experiment utilizing more bacterial clones from the same original evolution experiment (n = 222 clones, virulence relative to the average virulence of the clones evolved at 31 • C when incubated at 24 • C). For each model parameter, 95 % credible interval and median of the posterior are shown. clones, Figure 3C). This experiment also confirmed that the clones randomly chosen for sequencing were broadly representative of the larger pool of clones isolated from 142 the evolved populations. Decreased virulence of the bacteria in the insect host accompanied with decreased prophage induction when cultivated outside of the host suggests 144 that prophage induction closely relates to the bacteria virulence in the insect host.
This was supported by the correlation between average strain virulence in waxmoth 146 larvae and average PP4 induction rates (Spearman's ρ = 0.52, p = 0.004).

148
To explore the mechanistic basis of the evolutionary changes in prophage induction rate and in bacteria virulence and the connection between those two traits, we analysed the 150 genomic data of the 28 evolved strains used in our study and of the reference strain preexisting the evolution experiment. Sequencing was performed with the PacBio SMRT 152 platform and is described in detail in a companion study (Bruneaux et al., 2021).
We did not find any known virulence factors among the identified proteins of PP4, 154 but it should be noted that 12 of its 44 predicted proteins were unannotated (Figure S1). Some accessory proteins found in P2-related phages are known to contribute 156 to bacterial virulence, protection against other bacteriophages or SOS prophage induction (Christie and Calendar, 2016). More generally, prophage-encoded virulence 158 factors are considered one of the benefits brought by prophages to their bacterial host explaining the maintenance of prophages in bacterial genomes (Fortier and Sekulovic,160 2013; Koskella and Brockhurst, 2014). Another possible connection between prophage induction and bacterial virulence is via the release of endotoxins upon cell lysis that 162 could be a causative agent of S. marcescens virulence since its lysates are known to be cytotoxic on their own (Petersen and Tisa, 2012). Structural proteins were well con-164 served between PP4 and related sequences, as is common for P2-like phages (Nilsson and Haggård-Ljungquist, 2007), with the exception of tail fiber proteins ( Figure S1).

166
Tail fiber proteins are known to be involved in phage host specificity (Scholl et al.,  Figure 4: Alignment of the genomes from the 29 sequenced strains showing the variable genetic loci. Each circular track represents a sequenced genome, for which the evolutionary treatment is color-coded. Minor alleles for genetic variants are shown on the genome tracks in light grey (SNP) and dark grey (indels). Ticks outside the last genome track indicate non-synonymous variants (i.e. non-synonymous SNPs and indels resulting in a frame shift). The outer line represents coordinates along the genome and the locations of the five predicted prophages. Prophages PP4 and PP7 are the prophages shown in Figure S1.  Supplementary Table S2 for details. B) Visualization of the association between phenotypic values and major (M) and minor (m) alleles of variant a and of the "pooled" variants for galactokinase and glycosyltransferase. Colors correspond to the evolutionary treatment applied to each strain.
All in all, 52 variable genetic loci were identified among the sequenced clones but none was located inside prophage PP4 sequence (Figure 4 and Supplementary Ta-172 ble S2; Bruneaux et al. (2021)). It can thus be reasonably expected that mutations elsewhere in the genome or epigenetic modifications should be responsible for differ-174 ences in prophage induction rates and bacteria virulence. To identify the molecular basis of those phenotypic changes, we investigated the association between phenotypic 176 traits and genetic variants present in at least two strains using Wilcoxon rank sum tests adjusted for false-discovery rate (Benjamini and Hochberg, 1995). We also used 178 the methylation data obtained from PacBio SMRT sequencing to test for association between phenotypic traits and adenosine methylation in GATC motifs (which are rec- Several of the variable genetic loci associated with prophage induction and bacteria 184 virulence pointed towards a potential role of changes in the biofilm structure and in the outer structure of the cellular envelope in modulating phage particle production and 186 virulence in the insect host (e.g. genes involved in peptidoglycan and LPS biosynthesis,      The heatmap shows Spearman's ρ between methylated fractions of variable m6A epiloci (rows) and phenotypes (columns; Ph., phage release; Vir., bacteria virulence in waxmoth larvae). Overlapping or closest (≤ 500 bp) downstream genes were assigned to each m6A epiloci. Probable gene functions were assigned to each gene product based on a manual literature search. Visualization of the relationship between m6A methylated fractions and phenotypic trait values is shown for the heatmap cells marked with a letter. The m6A epiloci used here were the variable m6A epiloci with a methylation fraction range ≥ 0.2 across sequenced samples and an uncorrected p-value ≤ 0.005 for Spearman's ρ with at least one phenotypic trait.
The functional annotation of genes associated with phenotypic changes via adeno- Overall, a striking feature of the genetic and epigenetic changes observed in our 212 study was that none of the associated genes was directly connected with thermal selection pressure, which was the primary selective pressure in the evolution experiment,

Origin of the strains used in this study
The reference strain was Serratia marcescens ATCC 13880. The evolution experiment  Table S1). Since the candidate prophage PP4 was the only one for which induction was detected in our preliminary assays, we focused our annotation and comparative 276 analysis efforts on this prophage.
The PHASTER output for PP4 contained a set of contiguous CDS in the S. marces-278 cens genome which were related to phage functions. It also provided putative locations of attL/attR attachment sequences but those were located between prophage CDS in-280 stead of being at their periphery. We searched for better attL/attR candidates by looking for the longest repeated motif located 3000 bp upstream and downstream of 282 the proposed prophage CDS. We found better candidates for attL/attR which encompassed all the PP4 CDS: the corresponding motif was 20-bp long, compared to the 284 12-bp long motif found by PHASTER. We defined PP4 as the sequence encompassed by those attL/attR sequences.

286
Annotation of PP4 CDS was manually curated by merging the annotation obtained for the S. marcescens genome (as described above) and the annotation provided by 288 PHASTER. A last attempt to identify the CDS for which the "hypothetical protein" status remained at this stage was performed by blasting their predicted protein se-

METHODS
Evolution in a phage-bacteria-insect system quences against the nr database using NCBI blastx server and its default settings, but excluding the Serratia taxid (blast run on 2020-04-10).

292
Identification of phages and prophages related to PP4 To identify phages related to PP4, we retrieved all the phage genomes available from GenBank record. Finally, we also searched the genome of the Serratia strain we used in our study for other prophages similar to PP4, using the same local tblastx approach 316 as described for KONIH1. The dot plot showed that PP7 was related to PP4, and we METHODS Evolution in a phage-bacteria-insect system included PP7 in our final comparisons.
Once we extracted the nucleotide sequences for phages P88, SEN4, KSP20 (two fragments available: AB452988.1 and AB452989.1) and for prophages KONIH1/1-3 320 and PP7, we identified the matches between their predicted proteins and those from PP4 by running a final local tblastx search of each of those sequences against the PP4 322 sequence.
Quantification of prophage induction using qPCR Methods, but a brief description is given below.

332
We used seven specific primer pairs targetting each of the candidate prophage regions and one additional primer pair targetting a chromosomal, non-prophage-related 334 bacterial gene (Supplementary Table S3) to quantify the amount of prophage DNA copies relative to the amount of bacterial genome copies present in a culture using of those were then incubated at 24 • C while the other ten were incubated at 31 • C.
Larva survival was monitored at 1-3 hour intervals by checking for larva movements, 400 and time of death was recorded as the inspection time when a larva was found dead.
Additionally, for each incubation temperature, ten larvae were injected with sterile 402 medium and ten with sterile water as controls. This setup was replicated four times, resulting in a total of 80 infected larvae per strain (40 to incubation at 24 • C and 40 to 404 incubation at 31 • C). Some larvae from the first replication block were discarded due to a technical problem, leaving three replication blocks instead of four for some strains.

406
We analysed the larva survival data using a Cox proportional hazards model, where replication block, larval body mass, culture optical density, strain identity, incubation 408 temperature and the interaction between strain identity and incubation temperature were included as fixed effects. In this type of model, the hazard function describes Epiloci were then associated with annotated genes: a gene was assigned to an 442 epilocus if the adenosine base was located within the gene coding region, or less than 500 base pairs upstream of the initiation codon in order to cover potential regulatory 444 regions of the gene. Several gene set approaches were then tested to try to detect biological functions or pathways related to the epiloci associated with phenotypic traits. 446 We used gene-ontology enrichment tests as implemented in the TopGO R package and KEGG pathway analysis with Wilcoxon rank-sum statistics to compare gene sets, but

Association between genetic changes and phenotypic traits 638
The variable loci most clearly associated with both phage release and bacteria virulence in the insect host were the variant a and the pooled variants related to galactokinase 640 (pooled variants b, c, and h) and related to glycosyltransferase (pooled variants e, 27, 28, and 29 ) ( Figure 5). These genetic variants were located in or close to (≤ 500 642 bp) genes annotated as transcriptional regulators (molybdenum-dependent transcriptional regulator and transcriptional regulator RcsB involved in motility and capsule and 644 biofilm formation in E. coli) and enzymes involved in the cell wall and outer membrane structure and biofilm formation (peptidoglycan synthase, two glycosyltransferases and 646 a cellulose biosynthesis protein BcsG) (Supplementary Table S2). Those genes point towards a potential role for modifications of biofilm structure and of the outer struc-648 ture of the cellular envelope in modulating phage particle production and virulence in the insect host. In particular, the three independent mutations located in a single Finally, we also noted that haplotype a, comprising eleven associated genetic loci, 656 was shared by 5 out of the 8 strains evolved at 38 • C and by the reference strain, but by none of the other sequenced strains. This points to the probable existence of some 658 standing genetic variation at the onset of the experiment, which was then subjected to selection during the experimental evolution (Bruneaux et al., 2021).

Association between epigenetic changes and phenotypic traits
In addition to nucleotide sequences, the data we obtained from the PacBio SMRT for methylation and regulatory proteins specific to the same region (Casadesús andLow, 2006, 2013), and can thus be subject to selection.

674
Among GATC motifs which were not fully methylated in our dataset, no association was found between evolutionary treatments and methylated fractions (Bruneaux et al., 2021). However, we identified adenosines for which changes in methylation level were associated with phenotypic changes in the traits measured here (phage induction 678 and virulence in an insect host). For a given phenotypic trait, GATC loci exhibit-ing both positive and negative correlations between methylated fractions and the trait 680 values could be observed (Figure 6, heatmap panel). Manual curation of the function of the genes associated with GATC motifs related to phenotypic changes showed 682 that many of them were involved in (1) transcription regulation, (2) nutrient capture and transport into the cell, (3) excretion into the outer medium, (5) biofilm forma-684 tion, adherence or motility, and (6) cell envelope structure (including peptidoglycan and lipopolysaccharide biosynthesis) (Figure 6, gene functions panel). Many of those 686 functional categories have been shown to be critical for pathogen virulence in other bacterial species, in particular for nutrient capture in the challenging host medium 688 (Ren et al., 2018;Liu et al., 2017), for recognition of the host habitat via its nutrient signature (López-Garrido et al., 2015;Krypotou et al., 2019) and for biofilm formation, 690 adherence and motility which have a key role in colonization and successful invasion of the host tissues (Turner et al., 2009;Luo et al., 2017). The numerous candidate genes 692 involved in lipopolysaccharide biosynthesis also suggest that the O antigen, which can classically be involved both in cell recognition by phages and in bacterial virulence 694 in its host (Chart et al., 1989;Li and Wang, 2012), could act as a major player of evolutionary trade-offs between bacterial virulence and resistance to phage infection.

Culture conditions for the temperature assays
Frozen stocks had been stored at −80 • C in 40 % glycerol, with evolved clones stored 700 in 100-well plates (Bioscreen measurement plates), in randomized order and reference clone stored in microcentrifuge tubes. A preculture step in 400 µl of SPL 1 % at 31 • C 702 was performed by using a cryo-replicator to inoculate evolved clones into a new 100well plate and by inoculating the reference strain into wells of another plate. After 704 24 hours, five identical 100-well assay plates containing both the 28 evolved clones of interest and the reference clone were prepared by transferring 40 µl of each preculture 706 into 360 µl of fresh SPL 1 % (1 well per clone, i.e. 29 wells occupied per plate). For the first day of assay, one plate was incubated at 31 • C, two plates at 24 • C and two plates 708 at 38 • C. After 24 hours, clones within a given plate were transferred to 29 previously empty new wells in the same plate (40 µl culture into 360 µl fresh medium). For the 710 second day of assay, the plate from 31 • C was kept at 31 • C, one plate from 24 • C was kept at 24 • C and the other was transferred to 38 • C, and one plate from 38 • C was kept 712 at 38 • C while the other was transferred to 38 • C. After 24 hours, plates were taken for sample processing. Extra wells containing sterile SPL 1 % medium were used on the 714 assay plates to monitor potential contamination during plate handling (which was not observed). The whole experiment was performed twice, starting with the same frozen 716 stocks but with independent precultures.

718
At the end of the second day of assay, each of the 29 cultures in each of the 5 assay plates was processed in the following way: 50 µl of native culture sample was transferred 720 to a 96-well PCR plate, while the rest of the culture was placed into a microcentrifuge tube, centrifuged at 13 500 g for 5 min and 50 µl of supernatant was transferred in the 722 96-well plate, resulting in two paired samples per culture (native and supernatant). Samples from a given assay plate were placed into the same 96-well plate. A DNase 724 treatment was then performed to digest DNA fragments which were not protected inside a bacteria cell or a phage particle. 5 µl of DNase I at 0.1 mg ml −1 were added 726 to each sample, followed by an incubation at 37 • C for 30 min. DNA was then released from bacteria cells and potential phage particles by incubating the samples at 95 • C 728 for 30 min after having added 5 µl of EGTA (20 mm, pH 8) in order to hinder DNase I activity. The sample plates were then stored at −20 • C until DNA quantification by 730 qPCR runs.
Quantification of DNA target sequences was performed using prophage-specific 732 primer pairs and one bacterial-gene-specific primer pair (Supplementary Table S3). Preliminary experiments using the reference strain at 31 • C having showed no de-734 tectable extra-cellular DNA at least for prophages 2 and 5, six qPCR were runs per 96-well sample plate from this experiment using primers for prophages 1, 3, 4, 6, 7 and 736 for bacterial gene purA2. Runs were performed using CFX Real Time PCR detection system (Bio-Rad laboratories, USA). Amplifications were performed in a final volume of 10 µl, containing 5 µl of 2 x IQ SYBR Green Supermix (Bio-Rad), 0.5 µl of forward and reverse primers (6 µm each) and 1 µl of undiluted sample. Amplifications for each 740 primer pair were performed on separate qPCR plates, with in-plate calibration samples for each run. Calibration samples were prepared by serial dilution of a stock solution 742 of purified Serratia marcescens DNA of known concentration, and ranged in concentration from 10 6 to 1 genome copy per qPCR well, based on the predicted molecular To simplify its presentation here, we will first explain the modelling part related to 764 the estimation of induction rates for each culture well, based on the Cq values for the native and supernatant samples obtained from qPCR runs with bacterial and prophage 766 primers, before explaining the incorporation of assay and evolutionary treatment effects.

768
Let c bact,nat be the number of bacterial chromosome copies present in a native sample. The value of c bact,nat is determined from the qPCR run using the bacterial-genespecific purA2 primers. Let c pro,nat be the number of prophage DNA copies present in the native sample for e.g. prophage KSP20. The value of c pro,nat is determined from the qPCR run using the prophage-specific primers. Let's assume that this prophage is induced into phage particles at an activation rate a, such that the number of phage particles present in the native sample c phg,nat is related to the number of bacteria cells (i.e. the number of bacteria chromosome copies) by c phg,nat = a × c bact,nat . Since the prophage primers can target the prophage sequence both in the bacterial genome and in phage particles, we have: c pro,nat = c bact,nat + c phg,nat c pro,nat = c bact,nat + a × c bact,nat c pro,nat = (1 + a) × c bact,nat After centrifugation, we assume most bacteria cells have been pelleted and most phage particles (if any) have remained in suspension. Let k be the concentration factor during centrifugation for this culture, so that k = (c bact,sup /c bact,nat ) where c bact,sup is the number of bacterial chromosome copies present in the supernatant samples, as determined by qPCR with purA2 primers (0 ≤ k ≤ 1). If c pro,sup is the number of prophage DNA copies in the supernatant sample determined by qPCR with the prophage primers and c phg,sup is the number of phage particles in the supernatant sample, and if we assume c phg,sup = c phg,nat (i.e. we assume that the amount of phage particles pelleted during centrifugation is negligible), we have: Thus, to summarize, the two fundamental equations that relate the four qPCR measurements for a given culture (c bact,nat / c bact,sup / c pro,nat / c pro,sup ) and the prophage activation rate a in this culture are: c pro,nat = (1 + a) × c bact,nat c pro,sup = c bact,sup c bact,nat + a × c bact,nat We describe below the integrated Bayesian model used to estimate phage activation rates based on those two equations and on the Cq data obtained from qPCR runs.

770
Note that in the model description below, all parameters corresponding to DNA concentrations are expressed in number of target copies per qPCR well (copies/well).

772
The model relates Cq values to DNA concentrations, using plate-specific calibration parameters (calibration samples were present in all qPCR plates). Firstly, the model likelihood component due to the calibration samples is (with n cal being the total number of qPCR wells containing a calibration sample in our dataset): where cal i is the expected number of target copies in the well (between 1 and 10 5 in our experiment), c well i is the actual number of target copies in the well, α [.] and β [.] 774 are calibration parameters describing the relationship between Cq values and DNA concentrations, run i ∈ {1 . . . n runs } is the index of the qPCR plate corresponding to 776 the calibration sample and Cq cal i is the observed Cq value for the sample. Note that α [.] and β [.] are plate-specific (i.e. they are indexed by run i ) to account for plate variability 778 in the qPCR efficiency, while the parameter σ Cq which accounts for the experimental noise in the observed Cq values is shared across all qPCR runs. Note also that we 780 use a Poisson distribution for c well i to more accurately describe the sampling process happening when pipetting the template from their preparative tubes into the qPCR 782 wells, especially at low template concentrations.
Secondly, we describe the model likelihood component due to the qPCR wells containing the experimental samples of unknown concentrations prepared from the cultures in the assay plates. For this, we set (with n unkn being the number of qPCR wells with samples of unknown concentration and cult i the index of the original culture for each unknown sample): c bact,nat,cult i or c bact,sup,cult i or c pro,nat,cult i or c pro,sup,cult i depending on whether the unknown sample is run with purA2 (c bact,.,. ) or prophage (c pro,.,. ) primers and whether it is native (c .,nat,. ) or from supernatant (c .,sup,. ). The likelihood due to unknown samples is then of the same form as for the calibration samples: The remaining deterministic relationships of the model and the priors used for unknown parameters to estimate are: ∀i ∈ {1 . . . n runs }, α i ∼ normal(µ = 40, σ = 10) β i ∼ normal(µ = 3.5, σ = 2) σ Cq ∼ half-Cauchy(scale = 2.5) for the parameters of the qPCR calibration curve for each run (note that σ Cq is shared across all qPCR runs) and: ∀i ∈ {1 . . . n cultures }, log 10 (c bact,nat,i ) ∼ uniform(0, 20) log 10 (c bact,sup,i ) = log 10 (k i ) + log 10 (c bact,nat,i ) log 10 (c pro,nat,i ) = log 10 (1 + a i ) + log 10 (c bact,nat,i ) log 10 (c pro,sup,i ) = log 10 (k i + a i ) + log 10 (c bact,nat,i ) log 10 (k i ) ∼ half-Cauchy(scale = 2) log 10 (a i ) + 4 ∼ gamma(µ = 2, σ = 2) for the characteristics of a given culture well. Note that here, we assume that the 784 minimum value of activation rate a is 10 −4 , which is approximatively the lower sensitivity threshold predicted for our method when we assume that Cq values are measured 786 with a standard deviation σ Cq ≈ 0.48 (Supplementary Figure S7). We model this as (log 10 (a i ) + 4) following a Gamma distribution. In this explanation, we use fixed values 788 for the parameters of the Gamma distribution, but when we will introduce the effect of assay and evolutionary treatment the µ and σ parameters of this Gamma distribution 790 will depend on the treatments. This model formulation is sufficient to obtain posterior distributions for log 10 (a i ) for each culture well i in the assay plates. To model the effect of assay and evolutionary treatment, we extend the model by modifying the parameters of the previous prior for a i : log 10 (a i ) + 4 ∼ gamma(µ = 2, σ = 2) by: where assay i is the index of the assay treatment for culture i (assay i ∈ {1 . . . 5}) and str i is the index of the strain ID for culture i (str i ∈ {1 . . . 29}). The priors for the effect of assay treatments are: The strain effects include a hierarchical effect of the evolutionary treatment (four levels: three evolution environments plus the reference strain). The priors for the strain and evolutionary treatment effects are: (0, 10) where evo i is the index of the evolutionary treatment for strain i. [t j , t j+1 ) and is of duration dt j = t j+1 − t j . The survival data was transformed into a risk variable Y i (j) and an event count variable dN i (j) defined for all i ∈ {1 . . . N } and j ∈ {1 . . . T } by: The model assumes: where dλ 0 (j) is the increment in the integrated baseline hazard from t j to t j+1 and βz i is the product of the model parameters and of the covariate values for larva i. The term βz i corresponds to: where blk i , BM i , OD i , str i , and incub i are respectively the replication block, body mass, preculture OD, injected strain ID (str i ∈ {1 . . . 29}) and incubation temperature (0 for 24 • C, 1 for 31 • C) for larva i. Square brackets indicate indexing of a vector parameter; β blk is a vector containing the replication block effects and β str|incub24 and β str|incub31 are vectors containing the strain effects in the 24 • C and 31 • C incubations, respectively.
To model the effect of the evolutionary treatment, we set, for k ∈ {1 . . . 29}: where the vector evo allows to map the strain ID and one of the four evolutionary 802 treatments (three different temperature regimes plus reference strain). The priors we used were: (0, 10) and for all j ∈ {1 . . . T }: with c = 0.001 and r = 0.1. We used the first replication block and the effect of the reference strain in the 24 • C incubation as the references: We ran four chains in parallel with the JAGS MCMC sampler for 10 000 iterations 804 per chain, of which the first 5000 were discarded as burn-in. Model convergence and chain mixing was assessed by visual examination of trace plots and calculation ofR 806 values.
Selection of m6A in non-fully methylated GATC motifs 808 The method to identify GATC loci which were not fully methylated in our dataset was reported in a companion study (Bruneaux et al., 2021). Briefly, we calculated for each 810 GATC locus the distance between the point defined by its methylated fractions on the plus and minus strand and the point corresponding to full methylation on both strands

820
Supplementary Figure S2: Setup of the evolution experiment from which clones were isolated. One randomly selected clone per evolved population was used for sequencing. Open arrows after experimental evolution indicates steps where evolved clones were grown under common garden conditions. Details of the evolution experiment are available in Ketola et al. 2013(Ketola et al., 2013.
Supplementary Figure S4: Longevity of waxmoth larvae at two incubation temperatures after injection with experimental Serratia marcescens strains. Each vertical lane shows larvae injected with a given strain. The Medium and Water lanes show control larvae injected with sterile medium and sterile water, respectively. Longevity is corrected for the effect of replication blocks, culture optical density and larva body mass. Dots are individual larvae. Boxplots center line, median; box limits, upper and lower quartiles; whiskers extend to the most extreme data points which are no more than 1.5 times the interquartile range away from the box. . The shape of the trajectories depends on the activation rate of the prophage, i.e. on how many phage particles are present per bacteria cells in the native sample. The colored dots matching the colored predicted trajectories represent simulations of qPCR estimations which would be obtained as the centrifugation removes more and more bacteria cells from the supernatant, assuming a precision of the Cq values σ cq = 0.48 and triplicates qPCR measurements for each culture well, as was done in our experiment. As can be seen on the figure, the sensitivity threshold to detect phage particles decreases as the depletion of bacteria cells becomes more complete. However, even at k values of 10 −3 , activation rates of 10 −4 and lower are not distinguishable from the absence of induction. The red dots represent the results for a hypothetical culture, with the top-right dot representing the native sample and the bottom left dot representing the supernatant sample.