Environmentally triggered evolutionary cascade across trophic levels in an experimental phage-bacteria-insect system

2 Environmental changes can cause strong cascading effects in species communities due to altered biological interactions between species (Zarnetske et al., 4 2012). Highly specialized interactions arising from the co-evolution of hosts and parasites, such as bacteria and phages, and short generation times of these species 6 could rapidly lead to considerable evolutionary changes in their biotic interactions (Kerr, 2012; Buck and Ripple, 2017), with potential large-scale ramifications to 8 other trophic levels. Here we report experimental evidence of cascading environmental effects across trophic levels in an experimental system where phage10 bacteria coevolution in an abiotically altered environment cascaded on bacterial virulence in an insect host. We found that the lytic cycle of the temperate phage 12 KPS20 induced at low temperatures led to selection in the bacterial host Serratia marcescens that tempered the likelihood of triggering the phage’s lytic cycle. 14 These changes in S. marcescens concomitantly attenuated its virulence in an insect host, Galleria mellonella. In addition, our data suggests that this effect 16 is mediated by mutations and epigenetic modifications of bacterial genes moderating the onset of the temperate phage’s lytic cycle. Given the abundance of 18 temperate phages in bacterial genomes (Canchaya et al., 2003), the sensitivity of the onset of their lytic cycle to environmental conditions (Howard-Varona et al., 2

temperature over the last day of a given assay, rather than experiencing a temperature change between the two days: ending an assay at 31 • C induced about three times more 82 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. The consistency between higher 84 phage activation rates at lower assay temperature and the selection of lower phage activation rates in clones evolved at those temperatures suggests that the presence of 86 KPS20 prophage is having an effect on bacterial fitness, via phage release and cell lysis, especially at lower temperatures. We did not use 38 • C as the high incubation temperature since waxmoth larvae cannot 96 survive at this temperature. A Cox proportional hazards mixed model, controlling for the body mass of the larvae and the initial density of the bacterial sample (Figure 2A), 98 revealed that average virulence of clones evolved at hot temperature (38 • C) tended to be higher than for clones evolved at lower mean temperature when larvae were in-  Figure S4). This experiment also confirmed that the  Figure 2: Effect of evolutionary treatment on strains virulence in waxmoth larvae at two incubation temperatures. A) relative hazards for individual strains. The relative hazards were estimated from a Bayesian implementation of a Cox proportional-hazards model and are corrected for the effects of injection batch, larval body mass and optical density of injected cultures. All relative hazards are relative to the hazard rate from the stock strain in incubation at 24 • C (denoted by a broken horizontal line). B) Mean relative hazards per evolutionary treatment and per incubation temperature (exp(µ evo )) and standard deviation of the log−relative hazards per evolutionary treatment and per incubation temperature (σ evo ), as estimated by the model. For each variable, 95 % credible interval and median of the posterior are shown.
clones randomly chosen for sequencing were broadly representative of the larger pool of clones isolated from the evolved populations. 110 When put in relation with the phage activation results, decreased virulence in the insect host accompanied with decreased amount of virions when cultivated outside of 112 the host suggests that phage production by the bacteria closely relates to its virulence in the insect host. This was strongly supported by the correlation between average 114 strain virulence in waxmoth larvae and average KSP20 activation rates (Spearman's ρ = 0.52, p = 0.004). Using whole-genome sequencing, we identified 54 variable loci across all evolved 130 clones used in this study (n=28) compared to the stock strain (Figure 3 and Supplementary Tables S2 and S3). We investigated the association between genetic variants 132 present in at least two strains and phenotypic traits using t-tests and adjustment of The variable loci most associated with phage release were a, 29, 31 and 39 (fdrcorrected p-values between 0.01 and 0.1) and those most associated with virulence in 140 the insect host were a, 11, 28 and 29 (fdr-corrected p-values between 0.01 and 0.1) ( Figure 4). These genetic variants were located in or close to (< 500bp) genes anno-    Figure 5: Association between epigenetic changes and phenotypic traits. The heatmap depicts the largest correlations between m6A methylation fractions associated with a given gene (rows) and phenotypes (columns). Only m6A from GATC motifs which were not fully methylated across sequenced strains and which were significantly associated with at least one trait are shown. A white dot indicates a significant association (uncorrected p-value < 0.005 for Spearman's ρ). Probable gene functions are based on manual curation. A more detailed view of the genes associated with each trait is presented in Supplementary Tables S5-S11.   Table S3).
The strong trophic cascade caused by evolutionary adaptation of bacteria to its 212 phage, with consequences on virulence in insect host is a novel finding. As bacteria have a strong impact on the biosphere and biochemical cycles, and furthermore act as  Figure S2).

DNA extraction and sequencing 244
Selected clones (n = 10, 8 and 10 from 31 • C, 38 • C and fluctuating 24-38 • C evolutionary treatments, respectively) were grown from a frozen stock (40% glycerol) using a 246 cryoreplicator into 400 µL of SPL 1% overnight. A preculture of the stock strain was also initiated from a frozen aliquot in a similar volume. 150 ml of SPL 1% were inocu-248 lated with these precultures the next day and grown for 24 hours. Cells were pelleted and bacterial DNA was extracted using the Wizard Genomic DNA Purification Kit 250 from Promega (WI, USA), following the instructions provided by the manufacturer.
Extracted DNA was resuspended in water and one DNA sample (20 to 60 µg) per   Table S1).

Analysis of genetic variation 290
The chromosome consensus sequences of the 29 strains were aligned using Mugsy (Angiuoli and Salzberg, 2011). Variable loci were identified using a custom Python 292 script to identify variable positions in the alignment and to extract allelic information for each sequenced strain. To investigate the association between genetic variation 294 and phenotypic traits, we ran one t-test per genetic variant and phenotypic trait combination (using only genetic variants present in at least two strains). P-values were 296 corrected for multiple testing using the false-discovery rate method (Benjamini and Hochberg, 1995).

Analysis of epigenetic variation
Epigenetic data consisted of the methylation fraction for adenosine bases in all GATC 300 motifs present in the stock strain genome (38 150 GATC palindromes were present in the stock strain genome, corresponding to 76 300 adenosine bases for which methyla-  Epiloci were then associated with annotated genes: a gene was assigned to an epilocus if the adenosine base was located within the gene coding region, or less than 320 500 base pairs upstream of the initiation codon in order to cover potential regulatory regions of the gene. Several gene set approaches were then tested to try to detect bio-322 logical functions or pathways related to the epiloci associated with phenotypic traits. We used gene-ontology enrichment tests as implemented in the TopGO R package and 324 KEGG pathway analysis with Wilcoxon rank-sum statistics to compare gene sets, but mostly only very general biological functions were detected with those approaches, 326 such as amino acid or carbon metabolism, nutrient transport and translation (data not shown). Since those approaches are targetting the detection of changes affecting a given 328 biological function or pathway on average, but are not efficient to detect single genes which might affect phenotype, we decided to generate lists of top candidate genes as- were injected simultaneously; ten of those were then incubated at 24 • C while the other 382 ten were incubated at 31 • C. Larva survival was monitored at 1-3 hour intervals by checking for larva movements, and time of death was recorded as the inspection time 384 when a larva was found dead. Additionally, for each incubation temperature, ten larvae were injected with sterile medium and ten with sterile water as controls. This setup 386 METHODS Evolution in phage-bacteria-insect system was replicated four times, resulting in a total of 80 infected larvae per strain (40 to incubation at 24 • C and 40 to incubation at 31 • C). Due to some injection issues with 388 a clogged syringe in the first replication block, data from only three replication blocks were used for some strains. 390 We analysed the larva survival data using a Cox proportional hazards model, where replication block, larval body mass, culture optical density, strain identity, incubation  sample processing. Extra wells containing sterile SPL 1 % medium were used on the assay plates to monitor potential contamination during plate handling (which was not 578 observed). The whole experiment was performed twice, starting with the same frozen stocks but with independent precultures.

580
Sample processing and qPCR runs 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: 582 50 µl of native culture sample was transferred 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 584 and 50 µl of supernatant was transferred in the 96-well plate, resulting in two paired samples per culture (native and supernatant). Samples from a given assay plate were 586 placed into the same 96-well plate. A DNase treatment was then performed to digest DNA fragments which were not protected inside a bacteria cell or a phage particle. 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 (k ≤ 1), 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. 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 and that are used in the Bayesian model are: The priors and the deterministic relations of the statistical model 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 (k i ) ∼ half-Cauchy(scale = 2) log 10 (a i ) + 4 ∼ gamma(µ = 2, σ = 2) 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 ) for the characteristics of a given culture well. Note that here, we assume that the 630 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 632 with a standard deviation σ Cq ≈ 0.48 (Supplementary Figure S6). We model this as (log 10 (a i ) + 4) following a Gamma distribution. In this explanation, we use fixed values 634 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 636 will depend on the treatments.
The likelihood of the model due to calibration samples, for which the concentration values cal i are known and the Cq values Cq cal i are observed, is (with n cal the number of qPCR wells with a calibration sample): where run i is the index of the qPCR run. For qPCR wells containing samples of unknown concentration prepared from the culture wells in the assay plates, we set (with n unkn the number of qPCR wells with unknown samples 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: 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 (1 ≤ assay i ≤ 5) and str i is the index of the strain ID for culture i (1 ≤ str i ≤ 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 stock strain). The priors for the strain and evolutionary treatment effects are: where evo i is the index of the evolutionary treatment for strain i. into T = 20 intervals, so that the s i,i∈{1...N } values were homogeneously distributed across intervals (i.e. all intervals contained approximatively the same number of death 644 events). Intervals were defined by their boundaries t j,j∈{1...T +1} , such that interval j is [t j , t j+1 ) and is of duration dt j = t j+1 − t j .

646
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 the vector evo allows to map the strain ID and one of the four evolutionary treatments (three different temperature regimes plus stock strain).