One-step growth curves

In order to calculate the burst size of DMS3, we performed one-step growth curves (54). Overnight, stationary phase PA14 cultures were diluted 1:100 and grown in LB until they reached an OD of 0.5. Phage was added to these cultures at a 1:10 volume ratio, for a final MOI of 0.01, mixed well, and incubated at 37°C for 5 minutes. In order to calculate adsorption, a fraction of the sample volume was immediately spun down for 5 minutes at 30,000 rpm, and the supernatant was stored with 10% (vol/vol) chloroform to quantify the remaining free phage. To begin the one-step growth curve, the remainder of the samples were added to a pre-warmed medium at a 1:100 ratio to stop adsorption and grown on a roller drum at 37°C for the next 2.5 hours. Samples were taken every 10 minutes and mixed with 10% (vol/vol) of chloroform for later quantification of PFUs.

We calculated burst size with the following equation: the total phage produced (the difference between the average PFUs before and after the burst) was divided by the number of cells that were infected (estimated by the number of adsorbed phages multiplied by the estimated number of cells that proceeded through lysis). Our measurements of an 8% lysogenization frequency, estimated by spot-on-lawn assay, corresponded to measures in reference (55) (data not shown).

Here, β is burst size; p₂ and p₁ are the second and first plateaus, respectively; $c_i$ is the number of infected cells; and $l$ is the lysogenization frequency. The burst size of DMS3 is 41.8 ± 8.4 phages per lysed cell (Fig. S2). All PFUs were enumerated by spotting the phage-containing fraction on 0.5% double agar overlay plates.

Spontaneous induction measurements

Because multiple inputs could contribute to a higher raw PFU value in the stationary phase and because the stationary phase itself could be an inducing condition for some viruses (including phage Mu) (56), we chose to measure spontaneous induction separately in both exponential and stationary phase, which necessitated the normalization of PFU values with the CFU values. This was also a desirable method for controlling whether increased growth rate was responsible for increased PFU values.

Growth curves were started from overnight cultures of replicate purified isolates. We washed cultures three times by resuspension in fresh LB media, normalized the OD to 0.2, and diluted them 1,000-fold. Time points were taken at 0 and 2–7 hours to capture exponential phase growth and plated for both CFUs and PFUs. Samples were incubated at 37°C on a roller drum. We measured spontaneous induction ($q$) by taking the difference in the number of viral particles released by cells growing in the exponential phase and normalizing by the estimated burst size, the average growth of the culture, and the total time the culture grew.

Here, $q$ is spontaneous induction and has units of burst cells per time; $∆V$ is the total increase in virion particles; $\beta$ is the burst size; $∆t$ is the change in time; $C$ is the average amount of cells in the culture. To account for exponential growth, all calculations were based on the linear regression of the log₁₀-transformed $∆V$ and $C$ values.

We also calculated spontaneous induction based on calculations in reference (57). This paper approximated spontaneous induction at each time point with this formula: ${\displaystyle q=\frac{V}{\beta C}}$ ; here, $V$ is the number of virion particles at that time; $\beta$ is the burst size; and $C$ is the average amount of CFUs at that time. These methods produced qualitatively similar results (Fig. S3).

To measure induction in the stationary phase, we began growth curves in the same way as the exponential phase measurements. Time points were taken at 0, 10, 12, 15, and 18 hours and plated for both CFUs and PFUs. To calculate the rate of spontaneous induction, we used the same formula as above, except with un-regressed logged values.

Genome sequencing

In order to create a viral reference, DMS3 DNA was extracted from filtered Lys2 overnight supernatant using the Phage DNA Isolation Kit (Norgen Biotek Corp, Cat #: 46800) following the manufacturer’s instructions. Libraries were prepared with a Biomek 4000 liquid handler (Beckman-Coulter). We quantitated libraries with a Qubit fluorometer (Life Technologies Corporation, REF #: Q32866). Libraries were pooled and submitted for 2 × 250 paired-end sequencing by the Roy J. Carver Biotechnology Center at the University of Illinois Urbana-Champaign with an Illumina NovaSeq 6000. We received ~3.8 million reads with about 100× coverage.

We inoculated evolved isolates and ancestral controls in 2 mL deep well plates in LB and grew them overnight. We extracted gDNA with the Beckman-Coulter gDNA extraction kit as above using the Nextera Flex Library Preparation Kit (Illumina). We quantitated libraries with a Qubit fluorometer (Life Technologies Corporation, REF #: Q32866). Libraries were pooled and submitted for 2 × 250 paired-end sequencing by the Roy J. Carver Biotechnology Center at the University of Illinois Urbana-Champaign with an Illumina NovaSeq 6000. We received an average of about 5 million reads per genome. All raw reads are available on the NCBI database under BioProject number PRJNA1021667.

Genome analysis

We ran a custom QC pipeline on our raw FASTQ reads, available on Github (http://www.github.com/igoh-illinois). Briefly, the Illumina adaptor sequences were trimmed using TrimmomaticPE version 0. Read quality was checked with FastQC version 0.11.9 (options: --noextract -k 5 -f fastq). Reads were aligned using BWA-MEM (58) with default options to a two-contig reference genome containing both our reference PA14 sequence and our reference DMS3 sequence. To identify chromosomal mutations, we ran Breseq (59) with default settings on the trimmed and quality-controlled reads. SAM files were checked manually in both IGV version 2.12.3 (60) and Tablet version 1.21.02.08 (61).

To identify insertion sites of transposable phage, we ran a second, custom pipeline, available on Github (http://www.github.com/igoh-illinois). The pipeline identifies insertion sites at nucleotide resolution by identifying reads that map to both the host and viral chromosome (“split” reads). It further identifies insertion sites by finding reads that have been split on either side of a 5-bp window, creating a small overlapping region when mapped back to the host genome. This is the result of a 5-bp duplication, which is characteristic of Mu and Mu-like phage transposition (62, 63).

While manually verifying the phage insertion sites in IGV, we found putative large duplicated regions of the host chromosome. To verify these duplicated regions, we used the depth command in Samtools (64) to find the number of reads that covered each position in the genome and graphed this using RStudio (R version 4.3.1) (65). In order to assess the protein content of deleted and duplicated regions, FASTA files of sequence of the reference PA14 genome were given to the eggNOG-mapper-v2 pipeline (settings: genomic data; default options) (66).

We used the CRISPR Comparison Tool Kit (CCTK version 1.0.0) to identify and compare CRISPR arrays of the lysogens (settings: crisprdiff; default options) (67).

Mitomycin C induction experiments

Single colonies of the strains of interest were inoculated in LB, grown overnight at 37°C on a roller drum, and subcultured until they reached an OD of 0.5. Cultures were normalized, split, and incubated with or without 0.5 µg/µL mitomycin C for 3.5 hours, after which CFUs and PFUs were enumerated.

Model information

We use a compartmental model based on a system of ordinary differential equations. There are six lysogen compartments, each representing a lysogen characterized by a distinct rate of spontaneous induction, but otherwise being identical. Lysogens are induced at their associated rates and transition into the lytic state, which is followed by phage production and bursting. The model is along the lines of the ones presented and analyzed in our previous works (31, 68, 69). All lysogens are assumed to grow at the same rate, and the total bacterial population grows logistically. The model equations read:

We have partially non-dimensionalized the model so that 1 time unit in the simulations corresponds to about 50 minutes. The lysogeny growth rate is denoted by $r$, the spontaneous induction rates by $q_i,i=1,\dots ,6$, the rate of phage production by $\delta$, the burst size is by $\beta$, and the rate of viral degradation by $\mu$. The first two equations are decoupled from the last one describing the phage (since we do not consider superinfections in this model). Therefore, the dynamics of the bacterial compartments resemble those of generalized Lotka-Volterra competition.

Lysogen populations maintain diversity in CRISPR presence and function

Sequencing of lysogens from each population showed that lysogens evolved the CRISPR locus through a combination of mutation and large deletions and other structural variants (Table 1; Fig. 2). These mutations did not overlap with uninfected evolved populations, which exhibited point mutations in flagellar and quorum sensing loci (Fig. S1; Table S3), typical of other laboratory evolution experiments (72). Thirty-nine percentage (7/18) of isolates mutated or deleted spacer 1 in the CRISPR2 array, which contains the mismatched spacer that targets the integrated DMS3 (Fig. 2b; Fig. S5). Two mutations occurred in parallel between two different replicate populations in the evolved lysogen treatment: an A to G point mutation in the seed region of the self-targeting spacer, and an exact deletion of the self-targeting spacer and its upstream repeat (Fig. 2b; Fig. S6). Twenty-two percentage (4/18) of isolates had disruptions (three independent frameshift mutations and a small deletion) within the cas genes cas7 and cas8, which form part of the complex that mediates interference in P. aeruginosa (73). Of these three frameshift mutations, one led to the predicted loss of the cas7 RNA-binding domain, and two are predicted to interfere with the cas3 recruitment or interaction domain of cas8 (Fig. 2b; Table 1). One strain was recovered with a 3-kb deletion in the cas gene region that spans cas5 (also part of the interference complex) and cas7. While each strain with mutations and indels in the CRISPR array likely no longer targets the DMS3 prophage, they maintain CRISPR function. We confirmed that the evolved lysogens did not acquire new spacers (Fig. S6).

TABLE 1

Description of mutations recovered in evolved lysogens

Group	Sample ID	Mutated region	Description	Due to phage	Location on PA14-REF chromosome
ev_pop1	ev_pop1_1	CRISPR2 sp1-8	Deletion	No	2,936,222–2,936,676
ev_pop1	ev_pop1_1	tssK	Type 6 secretion system baseplate gene; nonsynonymous	No	94,370
ev_pop1	ev_pop1_2	CRISPR2 sp1-4	Deletion	No	2,936,462–2,936,675
ev_pop1	ev_pop1_3	243,737 bp	Includes CRISPR region	Yes	2,835,111–3,078,848
ev_pop1	ev_pop1_4	186,866 bp	Includes CRISPR region	Yes	2,839,227–3,026,093
ev_pop1	ev_pop1_4	27,769 bp	Duplicated region	Yes	1,120,309–1,148,078
ev_pop1	ev_pop1_4	cysT	Sulfate transport protein; (CAG)4→3	No	329,528
ev_pop1	ev_pop1_5	CRISPR2 sp1	A < G mutation in the second nucleotide	No	2,936,674
ev_pop1	ev_pop1_6	CRISPR2 sp1	GAT < G deletion of first and second nucleotide	No	2,936,673
ev_pop2	ev_pop2_1	88,123 bp	Includes CRISPR region	Yes	2,921,860
ev_pop2	ev_pop2_2	cas8	Frameshift (predicted loss of C-term helical bundle region)	No	2,929,846
ev_pop2	ev_pop2_3	194,512 bp	Includes CRISPR region	Yes	2,890,191–3,084,703
ev_pop2	ev_pop2_3	187,901 bp	Duplicated region	Yes	6,222,744–6,410,645
ev_pop2	ev_pop2_4	CRISPR2 sp1	Deletion	No	2,936,643–2,936,675
ev_pop2	ev_pop2_5	CRISPR2 sp1	A < G mutation in the second nucleotide	No	2,936,674
ev_pop2	ev_pop2_6	335,331 bp	Includes CRISPR region	Yes	2,811,042–3,146,373
ev_pop2	ev_pop2_6	244,019 bp	Duplicated region	Yes	876,284–1,120,303
ev_pop2	ev_pop2_6	Intergenic region	(GCCAAC)11→8	No	3,792,568
ev_pop3	ev_pop3_1	cas8	Frameshift (retention of first 57/435 aa)	No	2,930,725
ev_pop3	ev_pop3_2	CRISPR2 sp1	Deletion	No	2,936,644–2,936,675
ev_pop3	ev_pop3_3	cas7	Frameshift (predicted loss of RNA-binding domain)	No	2,927,776
ev_pop3	ev_pop3_4	60,538 bp	Begins in cas3 gene and includes CRISPR2 array	Yes	2,934,007–2,994,541
ev_pop3	ev_pop3_5	34,249 bp	Includes CRISPR region	Yes	2,903,549–2,937,794
ev_pop3	ev_pop3_6	Cas gene deletion	Deletion of cas7 and cas5, partial deletion of cas6 and cas8	No	2,927,288–2,930,093

Open in a separate window

Open in a separate window

Fig 2

Mutations in infected evolved strains are distinct from uninfected evolved strains. (a) Graph of mutations found in all infected evolved samples. (b) Close-up of point mutations and small deletions in infected evolved strains. The majority are found in the self-targeting spacer 1 in the second CRISPR array. Gray dashed lines indicate the location in the genes. (c) Close-up of large deletions in infected evolved strains. Gray dashed double lines indicate the boundaries of the CRISPR-Cas region. In panels a–c, the y-axis indicates sample ID; x-axis indicates the position on the PA14 reference chromosome. Points represent point mutations; triangles spanned by a segment represent deletions of the spanned region; squares spanned by a segment represent duplications of the spanned region. White fill indicates a mutation that was caused by a virus.

We observed that the remaining 44% (8/18) of evolved lysogen isolates carried deletions of varying sizes; the smallest being a 3-kb deletion in the cas genes, and the largest deletion, 335,331 bp, comprising about 5% of the PA14 chromosome. Seven deletions were centered on the CRISPR2 array (Fig. 2c). These deletions ranged from approximately 34 to 335 kb (mean = 163.3 ± 100.3 kb). Deletions were independent in each isolate, with no shared deletion boundaries even within the same culture (Fig. 2c). Similarly, a previous study showed that PA14, when challenged with free DMS3 virions and subjected to a short-term evolution experiment, evolved genome deletions encompassing the CRISPR region in the presence of self-targeting [see Extended Data Table 1 in reference (44)]. Therefore, in the evolved lysogen populations, we observed extensive coexistence between combinations of CRISPR spacer mutations and large entire deletions of CRISPR, demonstrating the importance of phage infection in determining distinct evolutionary trajectories in isolates in the same environment.

Evolved lysogens with large deletions are polylysogens

While confirming that the evolved lysogens had retained the phage at its original integration site (Materials and Methods), we found that the boundaries of the deletions of the CRISPR regions were composed of reads that mapped to both the PA14 and DMS3 chromosomes (hereafter referred to as “split” reads), indicating that the large deletions in these seven isolates resulted from a DMS3 transposition event that occurred from within the chromosome (Fig. 2c). Therefore, we consider these deletions to be phage mediated. The regions of the phage chromosome to which the split reads mapped and the orientation of phage reads at the boundaries of the deletions in two samples (1_4 and 2_3) suggest that more than one phage genome may be inserted in the gap (Fig. 3d through e; Table 2). Notably, we did not recover any mutations in the phage chromosome. This work clearly shows that the mechanism of deletion is through phage transposition, which may arise from failed or partial induction events within a lysogen.

Open in a separate window

Fig 3

Location and characterization of large duplicated regions in three evolved lysogen isolates. (a) Coverage plots of the host genome. Dashed lines indicate evidence of host-phage boundary. X-axis is the position on the PA14 chromosome. Where possible, the x-axis has been truncated to view all possible insertion sites. Y-axis is the depth of coverage at that nucleotide. Cyan line represents the putative duplicated region in between two viral insertion sites. Gray line is the region between the ancestral insertion sites. Regions of coverage in the deletion region in Lys2_ev_pop1_4 and Lys2_ev_pop2_6 are from domains in a TpsA1 and TpsB1 protein (annotated as a filamentous hemagglutinin protein), and in the 3′ end of the Lys2_ev_pop2_6 deletion, a domain in OprB (annotated as a porin). Both were identified via BLASTX. (b) Cartoon of resulting genome architecture of the evolved lysogens as represented in (a), from top to bottom. Dashed lines indicate a putative host-phage boundary. Gray arrows represent the phage at its ancestral insertion site; green arrows represent phage at new insertion sites. Blue arrows represent duplicate host sequences. The direction of the arrow indicates 5′ to 3′, with 3′ ending at the tip of the arrowhead. Dashed lines represent the junctions between phage and host. Small red arrows indicate host reads; numbers above the arrow indicate read count at that site that supports that junction. Ellipses and question marks between phage genomes represent uncertainty. Read files for each sample and at each location can be found in supplemental Data.

The large deletions, though centered around the CRISPR locus, had different deletion boundaries. To assess the shared gene content in these regions, we used eggNOG-mapper to query the functional protein content that was lost (Materials and Methods). In addition to CRISPR, the deleted regions were enriched with genes from COG category S (genes of unknown function), which included Type 6 secretion system-related genes tssF and tssG, and multiple pyoverdine system genes, which are often lost during lung colonization (6, 74, 75). Several pyoverdine biosynthesis and transport genes (fpvA, pvdE, pvdH, pvdI, pvdM, pvdM, and pvdO) were deleted in six of the seven isolates with large deletions. PvdS, a regulator of pyoverdine biosynthesis genes (76), and pvdR, a component of the pyoverdine efflux transporter (77), were also deleted in five of the seven isolates with large deletions. The co-localization of these virulence and defense gene cassettes with CRISPR may contribute to variation in these regions (78, 79).

Evolved polylysogens contain large duplications

Three isolates containing phage-mediated large deletions in CRISPR also had large duplications elsewhere in the genome (27, 188, and 244 kb; mean = 153 ± 113 kb). In these cases, these regions were not deleted but were doubled in coverage (Fig. 3a through c). The location of these large duplications exactly corresponded to additional insertion sites, which were recovered by our pipeline (Materials and Methods). Due to the orientation of the PA14/DMS3 split reads at these insertion sites, which faced away from each other rather than toward each other, and due to the fact that the split reads at the boundaries of the duplicated regions only represented about 50% of the total coverage, we interpret these regions to be large duplications with a phage genome in the middle, as opposed to two independent viral insertions (Fig. 3d through f). As the duplicated regions were not centered around a shared core, they were almost completely non-overlapping in their gene content. Only one gene was duplicated in two of the three isolates (2_3 and 2_6), which was betT, a choline transporter known to accumulate mutations in clinical isolates from CF patients (80,–82). Broadly, genes from category H (coenzyme metabolism) were represented in all three isolates, as well as from P (inorganic ion metabolism) and S (genes of unknown function), as in the large deletions. Several genes annotated as part of the major facilitator superfamily, a class of membrane-associated transporter proteins, were also duplicated in two of the three isolates, as well as genes from the moa family, which have recently been implicated in biofilm formation (83). The fact that deletions and duplications caused by transposition are found in the same genome suggests that at least a small number of cells induced phage transposition into multiple regions of the chromosome but did not lyse (84, 85).

We noticed that two of these duplications (in 1_4 and 2_6) independently evolved a shared boundary six nucleotides apart (at positions 1,120,309 and 1,120,303, respectively) in an intergenic region between the 3′ end of a hypothetical protein and the 3′ end of an AraC transcriptional regulator. Accordingly, we investigated whether these new insertion sites shared any sequence similarities. An analysis of all new lysogen insertion sites (e.g., the deletion and duplication boundaries) using the motif-finding software MEME Suite did not return any motifs, either using MEME (searching for a motif in a 15-bp region centered on the insertion site) or MEME-ChIP (searching for centrally enriched motifs 250 bp around the insertion site) (86, 87).

Neither the isolates containing large deletions nor the ones containing large duplications differed in their growth from other evolved strains that did not have large structural variation (Fig. S8A). This suggests that, under these conditions, the fitness costs to deletions, duplications, or carrying additional copies of the phage in the chromosome are smaller than the fitness gains by removing self-targeting. Although these phage-mediated deletions of CRISPR represent the addition of one to two phage genomes to the lysogen chromosome, the spontaneous induction rate of these isolates remains reduced relative to the ancestral PA14 lysogen strain (Lys2) (Fig. 1; Fig. S8B). Additionally, although viral output does not change with phage genome copy number after challenge with mitomycin C, cell survival is significantly increased with increased phage genomes (Fig. S8C), suggesting a possible mechanism of viral interference leading to cell survival, which may also contribute to a decreased spontaneous induction.

We thank Santiago Elena at I2SysBio and the University of Valencia for helping conceptualize and establish initial experimental evolution studies on lysogen fitness. We are grateful to George O’Toole for discussions as well as strains of bacteria and phages described herein, including Lys2 and DMS3. We gratefully acknowledge Dr. Alan Collins for helpful discussions during the early stages of this project. We thank Whitaker lab members for their helpful discussions. We thank Alvaro Hernandez and Chris Wright of the Roy J. Carver Biotechnology Center for sequencing expertise, as well as Jeff Haas of the School of Integrative Biology for crucial help with data storage and server access.

This work is funded by grants from the Cystic Fibrosis Foundation and the Allen Distinguished Investigator Award from the Paul G. Allen Foundation to R.J.W. and the National Science Foundation grant DMS-1815764 to Z.R. This research is also a contribution of the GEMS Biology Integration Institute, funded by the National Science Foundation DBI Integration Institutes Program, Award #2022049.

L.C.S. and R.J.W. conceptualized the study. L.C.S. performed formal analysis. L.C.S. and R.J.W. acquired funding. L.C.S., Z.R., and A.C.S. performed the investigation. L.C.S., Z.R., J.H.C., and R.J.W. designed the methodology. J.C.V. and R.J.W. provided resources. J.C.V. provided software. L.C.S. and Z.R. visualized the study. L.C.S. and R.J.W. wrote the original draft and reviewed and edited the manuscript.