Experimental evolution for niche breadth in bacteriophage T4 highlights the importance of structural genes.

Ecologists have long studied the evolution of niche breadth, including how variability in environments can drive the evolution of specialism and generalism. This concept is of particular interest in viruses, where niche breadth evolution may explain viral disease emergence, or underlie the potential for therapeutic measures like phage therapy. Despite the significance and potential applications of virus-host interactions, the genetic determinants of niche breadth evolution remain underexplored in many bacteriophages. In this study, we present the results of an evolution experiment with a model bacteriophage system, Escherichia virus T4, in several host environments: exposure to Escherichia coli C, exposure to E. coli K-12, and exposure to both E. coli C and E. coli K-12. This experimental framework allowed us to investigate the phenotypic and molecular manifestations of niche breadth evolution. First, we show that selection on different hosts led to measurable changes in phage productivity in all experimental populations. Second, whole-genome sequencing of experimental populations revealed signatures of selection. Finally, clear and consistent patterns emerged across the host environments, especially the presence of new mutations in phage structural genes-genes encoding proteins that provide morphological and biophysical integrity to a virus. A comparison of mutations found across functional gene categories revealed that structural genes acquired significantly more mutations than other categories. Our findings suggest that structural genes are central determinants in bacteriophage niche breadth.

INTRODUCTION

Experimental evolution has served an important role in the study of niche breadth in microbes (Kassen, 2002). Many of these studies have focused on virus systems and experimental exposure to different (or new) hosts, as the type, range, and availability of hosts present in the environment is an important source of selection pressure (Cooper & Scott, 2001; Crill, Wichman, & Bull, 2000; Duffy, Burch, & Turner, 2007; Enav, Kirzner, Lindell, Mandel‐Gutfreund, & Béjà, 2018; Kutnjak, Elena, & Ravnikar, 2017; Morley, Mendiola, & Turner, 2015; Novella et al., 1995). This and other experimental work in this realm has provided varying levels of support for the trade‐off hypothesis—that the ability to broadly exploit hosts (generalism) comes at the cost of the ability to exploit any single niche (specialism) (Bedhomme, Lafforgue, & Elena, 2012; Novella, Hershey, Escarmis, Domingo, & Holland, 1999; Turner & Elena, 2000). In particular, there have been relatively few rigorous treatments of the evolutionary genomics of niche breadth expansion, where the genetic determinants of niche breadth are resolved with functional inferences drawn between mutations and phenotypes. Such studies continue to have practical relevance for questions about the phenotypic manifestations of niche breadth evolution, whether trade‐offs arise as a consequence of adaptation, and about the molecular signatures of such evolution.

In this study, we sought to characterize the phenotypic and molecular changes associated with niche breadth evolution in a model bacteriophage. We employ Escherichia virus T4 as a system for experimental evolution. T4 is of particular interest because it is among the most well‐studied and fully characterized viruses. T4 is also surprisingly complex, with a genome ~170 kb in size and a gene density four times greater than that of herpes viruses and twice that of Escherichia coli. Of its 300 genes, approximately 160 have been functionally characterized (Miller et al., 2003), which provides the opportunity to resolve molecular mechanisms responsible for niche breadth evolution. We evolved T4 for 50 generations in one of three host environments: (a) E. coli C, (b) E. coli K‐12, and (c) daily alternation of E. coli C and E. coli K‐12 (Figure 1).

Figure 1

Experimental evolution schematic. 15 populations were seeded with a previously isolated T4 ancestral clone and split into three evolutionary histories with different host environments: Five populations were exposed to Escherichia coli C; another five were exposed to E. coli K‐12; and the last five were exposed to E. coli C and K‐12 in daily alternation. Serial passaging occurred for 20 days, which is equivalent to approximately 50 generations. Plaque assays on the original host E. coli B and the selection host E. coli C and/or K‐12 were performed for quantification of evolved phage samples on passages 6, 12, and 18. To measure productivity, assays on E. coli C, K‐12, and B were performed for the same evolved samples. Following the productivity assay, samples were quantified using E. coli B, which generated the final estimate of productivity (log10 titer [pfu/ml]). Sequencing was performed on the complete genomes of the T4 ancestor and 15 evolved populations at passage 18

Selection on either E. coli C or E. coli K‐12 mimics a constant environment, which is predicted to drive the evolution of specialists; whereas selection on the alternating hosts mimics a temporally variable environment, which is predicted to drive the evolution of generalists (Turner & Elena, 2000). Our results reveal the complexity of niche breadth evolution, with some populations demonstrating properties of a trade‐off, others less so. The genomic data reflected patterns across evolutionary histories: new mutations were overrepresented in genes that encode structural virion proteins. Notably, this pattern implies that structural genes—and in particular, those that function in host recognition, infection, and stability—are important in niche breadth evolution, regardless of conditions that promoted a particular ecological strategy (specialism or generalism). We discuss these findings in detail, and reflect on their implications for general viral ecology, and for the various arenas where bacteriophage niche breadth has practical utility—in disease emergence, public health surveillance, and efforts to engineer bacteriophage for therapeutic purposes.

MATERIALS AND METHODS

T4 and bacterial strains

This study used Escherichia virus T4 (American Type Culture Collection [ATCC] #11303‐B4) and three wild‐type bacterial hosts: E. coli B (ATTC #11303), E. coli C (Coli Genetic Stock Center 3,121), and E. coli K‐12 (Coli Genetic Stock Center 4,401). T4 infection of E. coli K‐12 has been well documented (Yu & Mizushima, 1982), but not E. coli C, which is a strain normally used for the propagation of Escherichia virus phiX174 (Wichman, Millstein, & Bull, 2005). E. coli B is the strain that has been historically used for the propagation of T4 (Demerec & Fano, 1944) and currently recommended by the ATCC for T4 propagation. Bacteria were stored as 25% glycerol stocks at −80°C; isolated bacterial colonies were obtained by streaking stocks onto lysogeny broth (LB) agar (Miller formula) petri dishes (BD, VWR International). All assays and serial passaging utilized liquid bacterial cultures which were prepared daily; single colonies were inoculated into glass culture tubes containing 5 ml of LB broth and incubated overnight at approximately 180 RPM (C1 Platform Shaker, New Brunswick Scientific) and 37°C. The T4 ancestral clone used to initiate experimentally evolved lines was isolated by plaque purification on LB agar petri dishes containing overlays comprised of 10 µl serially diluted phage lysate, 30 µl (~108 colony‐forming units [CFU]/ml) bacterial host culture, 270 µl LB broth, and 4 ml LB soft agar (7.5 g/L agar). Following overnight incubation at 37°C, a single well‐formed plaque was isolated, saturated with 500 µl LB broth, treated with 4% chloroform (VWR Life Science), and stored at 4°C.

Experimental evolution of T4

To initiate the experimental evolution, the T4 ancestral clone was used to seed 15 populations. Three sets of five parallel populations were evolved in one of three host environments: (a) E. coli C, (b) E. coli K‐12, and (c) a daily alternating pattern of E. coli C and E. coli K‐12 (Figure 1). The initial passage was perf(ormed in LB agar petri dishes containing overlays comprised of the T4 ancestral clone, appropriate bacterial host, 270 µl LB broth, and 4 ml LB soft agar, with initial infection occurring at MOI = 0.001 (volumes for T4 and bacteria were variable due to differences in host growth). Following overnight incubation at 37°C, phage lysates were collected, treated with 4% chloroform, centrifuged at 7,500 rpm and 4°C for 20 min (J6‐MI, Beckman Coulter), and sterilized using 0.2 µm syringe filters (Pall Life Sciences). 500 µl of each phage lysate was saved for storage at 4°C. These steps ensured that all bacteria were removed from phage lysates prior to usage for subsequent passages. Passages following the initial infection were plated as overlays on LB agar petri dishes, which contained 10 µl of phage lysate, 30 µl bacterial host culture, 270 µl LB broth, and 4 ml of LB soft agar, followed by overnight incubation at 37°C. Lysates generated over the course of the experimental evolution were collected, purified, and stored in the same manner as the initial lysate. In all passages, phage populations were propagated with hosts derived from newly prepared overnight bacterial cultures (no possibility for coevolution). To prevent cross‐contamination between evolutionary histories, each set of parallel populations was confined to separate laboratory rooms, and lysate supernatants were collected and purified using only materials (i.e., micropipettes and plastics) contained to the same room. Serial passaging was performed for 20 days, which is assumed to be equivalent to ~ 50 generations of T4 growth (we elaborate on generation time in the Appendix A).

Phage quantification

Quantification was accomplished using plaque assays in which 4 µl of 10‐fold serial dilutions of phage lysates were spotted on LB agar petri dishes containing overlays comprised of 30 µl bacterial host culture, 270 µl LB broth, and 4 ml LB soft agar. Petri dishes were incubated overnight at 37°C, with plaque visualization and counting occurring the following day. Quantified phage samples were expressed as PFU/ml; each plaque was assumed to have originated from a single infecting phage particle; thus, one plaque was equivalent to one PFU. Throughout the experiment, E. coli B was used for quantification because it is a consistent and highly permissive host. Following the experimental evolution, three replicate plaque assays (n = 3) were performed on the 15 populations at passage 6, 12, and 18, using both the selection host/hosts and the original host, to confirm that E. coli B remained a more sensitive host. This yielded 315 data points ([5 populations × 2 single host histories × 3 time points × 2 quantification hosts × 3 replications] + [5 populations × 1 alternating host history × 3 time points × 3 quantification hosts × 3 replications]). Replicate titers were log10‐transformed and compared for differences due to quantification host using unpaired t tests; each comparison was comprised of 30 measurements because titers were pooled by evolutionary history. In total, 12 comparisons were made using GraphPad Prism v. 7.0b (GraphPad Software, La Jolla California USA).

Productivity assays

Productivity assays were performed to measure phage titers produced on three hosts: E. coli C, E. coli K‐12, and E. coli B. Following the experimental evolution, assays were performed on the 15 populations at passage 6, 12, and 18. Productivity assays were standardized by infection of phage samples and bacteria at MOI = 0.001 and measured total number of progeny produced. Each assay was performed on LB agar petri dishes containing overlays comprised of 10 µl diluted phage lysate, 1 µl diluted bacterial host, 270 µl LB broth, and 4 ml LB soft agar. Following overnight incubation at 37°C, phage lysates were collected and purified using the same methods as those applied to the serial passages. Lastly, phage lysates were quantified for titer using plaque assays (as described above) and E. coli B as the quantification host.

In total, three replicate productivity assays (n = 3) were performed, which resulted in 405 titer measurements (5 populations × 3 evolutionary histories × 3 time points × 3 assay hosts × 3 replications). Titer measurements used in the linear mixed effects model were log10‐transformed to improve normality. The model was generated using the lmer function of the lme4 R package (Bates, Mächler, Bolker, & Walker, 2015). The initial model included replication as an additional random factor; however, this parameter was dropped because the variance estimate was extremely small (variance = 0.00142) when compared with the variance of the residual error (variance = 0.11195). To confirm goodness of fit and ensure that model assumptions were met, model residuals were visually inspected (Zuur, Ieno, & Elphick, 2010). Plotted residuals were normally distributed and revealed no concerning patterns of heterogeneity in variances. The analysis of variance table was calculated using the ANOVA function of the car R package, which performs Wald chi‐square tests for linear mixed effects models (Fox & Weisberg, 2011). Follow‐up analysis was performed using the contrast function of the lsmeans R package to obtain pairwise comparisons among LS means (Lenth, 2016). P values generated from this analysis were adjusted using the Holm method, which is designed to give strong control of the family‐wise error rate while retaining more power than the Bonferroni correction. The results of this analysis are available in Appendix B.

Additional productivity assays (n = 3) were performed for the T4 ancestor in order to gauge baseline titer on the three assay hosts. These data were excluded from the linear model and analyzed using one‐way ANOVA followed by Tukey's HSD test. Titer measurements of the ancestor were necessarily excluded from the linear model because by definition, the ancestor did not have an evolutionary history nor was it serially passaged. All above analyses were performed in R v. 3.3.3 (R Development Core Team, 2017).

Phage DNA extraction

Phage lysates were newly generated from the original evolved phage samples for library preparation (grown on E. coli B, a highly permissive host). This step was performed for two reasons: (a) to preserve the original 500 µl generated directly from the serial passaging, and 2) phage DNA extraction protocols generally require at least 1 ml of high titer (>109) phage lysate to generate high yield, clean DNA. Lysates were created from evolved samples at passage 18, the final time point that was assayed for productivity (Figure 1). These were plated as overlays on LB agar petri dishes, which contained 10 µl of evolved phage sample, 30 µl bacterial host culture, 270 µl LB broth, and 4 ml of LB soft agar. Following overnight incubation at 37°C, lysates were collected and purified in the same manner as with the serial passaging.

1 ml of each lysate was first treated with 12.5 µl 1M MgCl2 (VWR Life Science), 0.4 µl DNase 1 (2000 U/ml) (New England Biolabs), and 10 µl RNase A (10 mg/ml) (Thermo Fisher Scientific). Following brief vortexing and room temperature (RT) incubation for 30 min, the following reagents were added: 40 µl of 0.5 M EDTA (Corning), 2.5 µl Proteinase K (20 mg/ml) (Thermo Fisher Scientific), and 25 µl 10% SDS (VWR Life Science). The mixture was vortexed vigorously and incubated at 55°C for 60 min, with additional mixing occurring twice at 20 min intervals. Phage DNA was extracted using an equal amount of phenol:chloroform:isoamyl alcohol (25:24:1) (VWR Life Science) to lysate. Following inversion mixing and centrifugation at RT for 5 min at 13K, the aqueous layer was removed, and two further extractions were performed in the same manner. Phage DNA was precipitated by addition of 1 ml 90% ethanol (Pharmco) and 50 µl 3M sodium acetate (Corning). Following incubation on ice for 5 min, DNA was gently mixed and centrifuged at room temperature for 10 min at 13K. The DNA pellet was washed with 500 µl 70% ethanol and centrifuged at RT for 10 min at 13K. The ethanol was decanted, and the pellet was air‐dried for ~ 20 min. DNA was reconstituted in 50 µl TE buffer (Invitrogen) and checked for purity and concentration with Nanodrop 2000 (Thermo Scientific).

Library preparation and sequencing

Library preparation was executed according to previously published protocols (Baym et al., 2015). Minor protocol modifications were made and are thus noted. In module 1, DNA was standardized with Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific). In module 4, DNA size selection was performed twice (rather than once). First, 0.55 × magnetic beads were added to genomic DNA to remove large (> 600 bp) fragments. Following incubation at RT for 5 min, the tubes were placed on a magnetic stand, and the supernatant was removed and transferred to new PCR strip tubes. In the second size selection, 0.3 × magnetic beads were added to the supernatant to bind DNA fragments of desired length (~250–600 bp). Following incubation at RT for 5 min, tubes were placed on a magnetic stand, and the supernatant was removed and discarded. The beads were washed once with 200 µl 80% ethanol, allowed to air‐dry for ~ 20 min, and resuspended in 30 µl TE buffer. In module 4, DNA was quantified using Qubit dsDNA HS Assay Kit; libraries were pooled at 5 ng per sample, and fragment size was assessed with Agilent 4,200 TapeStation System. The final pooled library sample was sent to Bauer Core Facility (Harvard University, Cambridge, MA) for QPCR quality control and sequencing using Illumina NextSeq 500 Mid Cycle (2  150 bp).

Analysis of sequence data

Demultiplexed reads were trimmed for Nextera adapter sequences using Trimmomatic with default settings (Bolger, Lohse, & Usadel, 2014). The reads were aligned to the T4 reference genome (RefSeq accession no. GCF_000836945.1) using Breseq v. 0.31.1 (Deatherage & Barrick, 2014). For populations, the polymorphism mode was enabled, and no junction predictions were made. Typical coverage depth was over 1,000 × for most positions in the genome with over 95% of the reads mapping. The ancestral clone sample was replicated in the library preparation stage to account for sequencing errors and potential variation between the reference genome and the ancestral clone. Following sequencing, mutations were called separately for the two ancestral clone samples; these replicates revealed complete agreement between the called mutations. The ancestral mutations were subsequently filtered from all evolved populations in order to account for mutations that were present in the T4 ancestor but absent in the reference genome. The default Breseq configuration requires a mutation to be found in at least 5% of the reads; hence, only SNPs present in the populations at 5% or greater were considered. All polymorphisms and their frequencies are available in Appendix B.

In the genomic analysis, sequencing error bias presented barriers for the determination of polymorphisms. To remedy this, mutations with consistent statistical evidence of strand bias were removed in circumstances where the mutation would not be called when considering information from one strand versus the other. Specifically, mutations were filtered due to strand bias in cases where there existed at least 5% difference between the forward and reverse strand (preventing cases where a mutation is found in 4.9% on one strand, and 5.1% on the other one, but filtering out cases where a mutation is found in 1% in one strand, but 6% on the other). xNPs (multiple adjacent nucleotide polymorphisms) similarly presented challenges to the NGS analysis pipeline, as it did not recognize sequential mutations as one mutational event, but rather multiple events. This issue manifests in the context of calling amino acid changes, that is, when two adjacent mutations occur within one codon, two different amino acid changes are detected for each mutation, rather than one amino acid change which incorporates both nucleotide substitutions. In these cases, reads were visualized to ensure all putative xNPs occurred within the same read (as opposed to a situation where mutations occurred separately in different reads) to validate the presence of xNPs. The amino acid calls were subsequently manually modified to reflect the true amino acid change. A similar issue occurred with indels; each indel was reported as multiple mutational events rather than one. These situations were resolved in the same manner as xNPs. A final issue occurred in which a fixed ancestral mutation was detected as a new polymorphism, when it decreased in frequency due to a new mutation. This new mutation was not detected because the nucleotide substitution was coincidentally the same base as that of the reference sequence used for alignment. This error was confirmed and manually corrected by visualizing reads containing the relevant position. A list of corrected positions is available in Appendix B.

Analysis of mutation distribution

The number of nucleotide sites occupied by each functional category was determined using previous gene designations (Miller et al., 2003) and the T4 genomic sequence associated with that study (Genbank accession no. AF158101.6). From this, the mutation rate (number of mutations per nucleotide site) of each functional gene category was calculated for the 15 evolved populations. All polymorphisms detected at or above 5% were considered for this analysis. A SRH nonparametric two‐way ANOVA (Sokal & Rohlf, 2011) was performed to test the influence of factors evolutionary history, functional gene category, and their interactions on rate of mutation. This analysis was first performed for all mutations and then repeated separately for synonymous and nonsynonymous mutations. The analyses were executed using the Scheirer–Ray–Hare function of the rcompanion R package (Mangiafico, 2018). Follow‐up analysis of pairwise comparisons was performed using Dunn's test, which was executed using the dunnTest function of the FSA R package (Ogle, 2018). P values generated from the pairwise comparisons were adjusted using the Holm method. All above analyses were performed in R v. 3.3.3 (R Development Core Team, 2017).

RESULTS

The baseline productivity of the T4 ancestral clone differs according to host strain

Before discussing the results corresponding to the state of phage populations postexperimental evolution, we first examined the baseline productivity of bacteriophage T4 across the various hosts (E. coli C, E. coli K‐12, and E. coli B) used in the experimental evolution. The purpose of this analysis was to address questions about how ancestral populations performed across backgrounds, prior to the evolution experiment. Phage productivity was chosen as an estimate for gauging performance on a particular host because it represents the number of infectious progeny produced (per infected bacterium or per parent virion) from an infection (Abedon, Herschler, & Stopar, 2001; Kerr, Neuhauser, Bohannan, & Dean, 2006). The ancestral clone displayed differing levels of productivity on each host strain, which was expressed as mean titer (log10 plaque‐forming units [PFU]/ml). Productivity was predictably highest on the ancestral host E. coli B (10.37 log10 PFU with 0.14 standard deviation [SD]), slightly lower on E. coli K‐12 (9.77 log10 PFU, 0.06 SD), and lowest on E. coli C (6.10 log10 PFU, 0.20 SD) (Figure 2a). Analysis of these data using one‐way ANOVA followed by Tukey's Honest Significant Differences (HSD) test indicated that productivity differences were significant (p < .01). Notably, T4 productivity was slightly different but high on E. coli K‐12 and B, while productivity on E. coli C was several orders of magnitude (on the order of 104) lower than the others.

Figure 2

Productivity of T4 clones. (a) Ancestral clones (before evolution). Productivity of the T4 ancestral clone was measured on the experimental evolution hosts, Escherichia coli C and K‐12, and on the original host, E. coli B. Each bar represents mean productivity (log10 titer [pfu/ml]) measured over three replicate assays, and error bars indicate 95% confidence limits. Analysis using one‐way ANOVA followed by Tukey's HSD test indicated that all comparisons of mean titer were significant (p < .01). (b) Evolved clones. The 15 populations were quantified for phage titer on the original host E. coli B and the selection host(s) E. coli C and/or K‐12 at passages 6, 12, 18. Plaque assays for each population were performed with threefold replication. Panels indicate one of three evolutionary histories, and each point represents grand mean titer (log10 PFU/ml) of the five populations within each evolutionary history; error bars indicate 95% confidence limits of the grand mean. A series of t tests comparing titers pooled by evolutionary history indicated that all comparisons involving E. coli C and B yielded significant differences (with the higher titer estimate occurring on E. coli B), and all comparisons involving E. coli K‐12 and B yielded no significant differences (see main text and Table 1)

Table 1

p Valuesa of t tests comparing log10 titer of T4 populations after plaque quantification on the evolved host(s) versus Escherichia coli B

Evolved host	Comparison	Passage
		6	12	18
E. coli C	C vs. B	<0.0001	<0.0001	<0.0001
E. coli K−12	K−12 vs. B	0.6603	0.3759	0.0424
Alternating	C vs. B	<0.0001	<0.0001	<0.0001
	K−12 vs. B	0.2034	0.6674	0.1273

The Bonferroni correction was applied to the critical value (α = 0.004) to adjust for multiple comparisons.

After evolution, E. coli B remained the most sensitive phage quantification host

Following experimental evolution, samples required quantification of phage titer prior to being assayed for productivity. E. coli B is commonly used for T4 quantification because it is a consistent and highly permissive host. We anticipated a potential problem because selection on the other hosts (E. coli C and/or E. coli K‐12) could lead to correlated antagonistic changes (i.e., reduced infectivity) on the original host.

Plaque assays were performed on the 15 evolved populations at three time points (passage 6, 12, and 18) to investigate whether E. coli B remained a more sensitive quantification host than the selected host(s) (E. coli C and/or E. coli K‐12). Results from these plaque assays indicated that phage titer measurements, expressed as log10 PFU/mL, on E. coli B were either equivalent or greater than those on the selected host(s) (Figure 2b). Specifically, the grand mean titer for E. coli C evolved populations was higher on E. coli B than that on E. coli C for all time points, indicating that E. coli B was a more sensitive quantification host. The grand mean titer for E. coli K‐12 evolved populations was equivalent whether quantified on E. coli B or K‐12 for all time points, indicating that these two strains were equally sensitive as quantification hosts. The grand mean titer for evolved populations with alternating hosts recapitulated the results of the single host evolved populations. For all time points, E. coli B provided equivalent titer measurements when compared with E. coli K‐12 and more sensitive measurements when compared with E. coli C. Furthermore, a series of t tests conducted on titer measurements pooled by evolutionary history indicated that all comparisons involving E. coli C and E. coli B as quantification hosts yielded significant differences (with the higher titer estimate occurring on E. coli B) whereas all comparisons involving E. coli K‐12 and E. coli B yielded no significant differences (Table 1). These results confirm that, despite selection on different hosts, E. coli B remained the most highly permissive host and thus was used for quantification of phage titer throughout this study.

Niche breadth evolution resulted in changes in phage productivity over time

After 20 days of host selection, the 15 evolved populations (5  3 replicates) were assayed for productivity on three hosts (E. coli C, K‐12, and B) at three equidistant time points (passage 6, 12, and 18). We measured phage productivity at three time points and across histories and assay hosts, relative to the productivity of the ancestor. This experiment yielded a total of 405 productivity measurements (including n = 3 replicate assays). Mixed effects linear modeling fit by restricted maximum likelihood was used to test for the effects of experimental treatments—evolutionary history, assay host, and number of passages—on phage productivity. Evolutionary history, assay host, number of passages, and their interactions were specified as fixed factors, and population nested within history was specified as a random factor. We used the model to generate predicted least squared (LS) means for productivity at each level of evolutionary history, assay host, and number of passages (Figure 3). Due to the presence of interactions, a type III analysis of deviance table was calculated to determine the significance of each fixed parameter and their interactions in explaining variation in phage productivity (Table 2). This analysis revealed that all fixed factors and their interactions had a significant effect on productivity.

Figure 3

Model predictions of mean productivity. Predicted LS means for productivity at each level of evolutionary history, assay host, and number of passages (estimates were not provided for individual populations because population was specified as a random factor). Horizontal rows denote the assay host, and vertical columns and colors denote evolutionary history. Each point represents LS mean productivity (log10 titer [pfu/ml]), and error bars indicate 95% confidence limits. Dashed black lines indicate the productivity of the T4 ancestor on a particular assay host to provide a frame of reference

Table 2

Analysis of deviance tablea, a for a mixed effects linear model testing the effect of each fixed parameter and their interactions on productivity (log10 titer) of evolved T4 populations

Fixed effects	Chi‐sq	df	Pr(>chi‐sq)
Evolutionary history	14.898	2	0.0005
Assay host	1,004.209	2	<0.0001
Number of passages	45.812	2	<0.0001
History × assay	45.366	4	<0.0001
History × passages	28.224	4	<0.0001
Assay × passages	115.056	4	<0.0001
History × assay × passages	74.234	8	<0.0001

Type III Wald chi‐sq test.

In addition, we performed pairwise comparisons on predicted LS means for productivity between evolutionary histories for each level of assay host and number of passages (see Figure 3 for a graphical representation and Appendix B for numerical estimates). Focusing on the final time point (passage 18) and assay host E. coli C, all contrasts of predicted means for phage productivity between evolutionary histories were significant (p ≤ .001), indicating that by passage 18, evolutionary history led to different outcomes of productivity on E. coli C. Predicted productivity was highest for populations evolved on E. coli C (LS mean = 7.14, 95% CI = [6.93; 7.35]), lower for populations evolved on alternating hosts (LS mean = 6.63, 95% CI = [6.41; 6.84]), and lowest for populations evolved on E. coli K‐12 (LS mean = 5.76, 95% CI = [5.55; 5.97]). At passage 18 and assay host E. coli K‐12, all contrasts of predicted means for phage productivity between evolutionary histories were not significant (p > .1), indicating that evolutionary history had the same effect on productivity on E. coli K‐12 at this time point. Furthermore, all predicted productivities on E. coli K‐12 were lower than that of the ancestor (E. coli C LS mean = 8.93, 95% CI = [8.72; 9.14]; E. coli K‐12 LS mean = 9.00, 95% CI = [8.79; 9.21]; alternating host LS mean = 8.69, 95% CI = [8.48; 8.90]). At passage 18 and assay host E. coli B, contrasts of predicted means for phage productivity between the alternating host history and the E. coli C and E. coli K‐12 histories were significant (p < .01) while the difference between the predicted means of the two single host histories was not significant (p = .60). Interestingly, all predicted productivities on E. coli B were also lower than that of the ancestor, with the lowest productivity belonging to the alternating history (E. coli C LS mean = 9.14, 95% CI = [8.92; 9.35]; E. coli K‐12 LS mean = 9.21, 95% CI = [9.00; 9.43]; alternating host LS mean = 8.68, 95% CI = [8.47; 8.89]).

Genetic signatures are suggestive of positive selection in evolved populations

Following the experimental evolution, the genomes of the 15 evolved populations at passage 18 were sequenced to analyze mutations and their frequencies. The small genome size of T4 coupled with current sequencing technology yielded over 1,000 × coverage depth for most positions in the genome, with over 95% of the reads mapping unambiguously to the T4 reference genome; this enabled the reliable detection of polymorphisms occurring at 5% frequency or greater. A total of 174 mutations occurred across 145 loci, with a range of 3–30 mutations per population and 77.6% occurring as unique to a particular population (for a graphical representation see Figure 4, and for a complete list see Appendix B). Of the total number of mutations, 2.8% occurred in intergenic regions; 6.9% were indels; 82.2% were single nucleotide polymorphisms (SNPs); and 8.1% were multi‐nucleotide polymorphisms (xNPs). All xNPs were nonsynonymous, and among them, 14.7% were synonymous, and 85.3% were nonsynonymous. Of the total number of detected SNPs, we observed a dN/dS ratio of 5.8. If observed indels and xNPs are also classified as nonsynonymous (each instance is conservatively counted as one nonsynonymous mutation despite causing more than one amino acid change in most cases), this ratio increases to 7.0 (Figure 5 and Table 3). Across evolutionary histories, we observed that few mutations reached fixation, with the preponderance of those being nonsynonymous mutations. These signatures indicate that positive selection shaped phenotypic changes during experimental evolution.

Figure 4

The distribution of new mutations detected in the evolved phage populations. Panels indicate one of three evolutionary histories, and lines within panels represent particular populations within an evolutionary history; numbers for each horizontal line correspond to the numbered populations that underwent passaging. Each point represents a new mutation (relative to the ancestor) detected at ≥5% frequency. Point shape indicates type of mutation (filled triangles: nonsynonymous; empty triangles: synonymous; circles: indel; squares: xNPs (multiple adjacent nucleotide polymorphisms); and diamonds: intergenic). Colors indicate functional gene category; labels above the top panel and dashed lines denote relative positions of virion structural genes

Figure 5

Classification of mutations. Panels indicate one of three histories; categories within panels indicate the total number of observed (>5%) and fixed (>80%) mutations for the five parallel populations within each history. Bar colors represent the type of mutation for each category: intergenic, synonymous, and nonsynonymous

Table 3

The total number of observed and fixeda mutations in the 15 sequenced populations, grouped by evolutionary history and sorted by type. The percentage of mutations observed and fixed in a given type are shown in parentheses

History	Category	All	Nonsyn	Syn	Intergenic
Escherichia coli C	All	85	74 (87%)	9 (11%)	2 (2%)
	Fixed	12	10 (83%)	2 (17%)	0
Escherichia coli K−12	All	45	35 (78%)	7 (15%)	3 (7%)
	Fixed	1	1 (100%)	0	0
Alternating	All	44	39 (89%)	5 (11%)	0
	Fixed	8	8 (100%)	0	0

Mutations are classified as fixed if they reach a frequency >80%.

Mutations in structural genes were predominant in evolved T4 populations

Mutations during evolution tended to cluster in regions of the T4 genome containing genes that code for structural proteins. Furthermore, of the 62 polymorphisms that reached a frequency of 50% or greater, all but four occurred in structural genes. To formally test whether different regions of the genome experienced different rates of mutation, and if the rates also depended on the evolutionary history of the population, a Scheirer–Ray–Hare (SRH) nonparametric two‐way ANOVA (Sokal & Rohlf, 2011) was performed. First, the number of nucleotide sites that each functional category occupies in the genome was determined. From this, the mutation rate (number of mutations per nucleotide site) of each category was calculated for the 15 evolved populations. Evolutionary history, functional category, and their interaction were specified as factors in the SRH nonparametric two‐way ANOVA. This analysis showed that the rate of mutation was different across functional gene categories (p < .001). There was neither a significant effect for evolutionary history (p = .588) nor for interactions between functional category and history (p = .789).

To examine pairwise comparisons between functional categories with regards to mutation rate, we used a Dunn's test, which revealed that only comparisons involving the structural gene category were significant (p < .01). Specifically, this indicated that the mutation rate of the structural gene category was significantly higher (i.e., mutations were overrepresented in structural genes) than all other categories.

DISCUSSION

In this study, we performed experimental evolution using Escherichia virus T4 to examine the phenotypic and molecular manifestations of niche breadth evolution. Our findings demonstrate that niche breadth evolution led to measurable changes in phage productivity. For populations grown on E. coli C, productivity on E. coli C increased while productivity decreased on both E. coli B and K‐12, indicating that adaptation on E. coli C led to reduced performance on both the original and unselected hosts. These results conform to expectations, set by other seminal studies, where evolved populations improve performance in the selective environment at the expense of diminished performance in alternate environments (Turner & Elena, 2000). For the alternating host evolved populations, the magnitude of increased productivity on E. coli C was smaller than that of the populations that only experienced the same host, which suggests that exposure to two hosts and temporal variability limited adaptation on E. coli C. This result also conforms to theoretical expectations; specialists are predicted to evolve faster than generalists because the time to fixation of favorable alleles is shorter for specialists (Whitlock, 1996).

Some experimental populations deviated from the expectations of simple trade‐off models

In contrary to the results of the E. coli C‐exposed populations, the observed changes in productivity (on E. coli K‐12) for the E. coli K‐12 and alternating host populations offered surprising results. For example, these treatment populations demonstrated decreased productivity on E. coli K‐12, despite being experimentally exposed to this host. These results were not only unforeseen, but also deviate from theoretical expectations. We predicted that evolution on E. coli K‐12 would lead to higher productivity with the single host exposed populations performing better than the alternating host populations. Instead, we observed decreased E. coli K‐12 productivity with no significant differences between the single (E. coli K‐12) and alternating host (E. coli K‐12/ E. coli B) histories. Though empirical evidence from previous viral experimental evolution studies have contradicted theory and found no cost associated with generalism (Bedhomme et al., 2012; Novella et al., 1999; Turner & Elena, 2000), our results offer a unique departure from past findings. This raises an important question specific to our system: why would T4 passaged on E. coli K‐12 evolve lower productivity on its selective host? One simple explanation is that productivity did not improve because it was not the target of selection, which undermines the neat (and convenient) equivalence between productivity and reproductive fitness. In contexts where competition is more synonymous with fitness, one might expect productivity to be completely decoupled from evolution. Indeed, prior studies have suggested that productivity can be negatively correlated with competitive ability in T4 in the setting of experimental evolution (Kerr et al., 2006).

Study limitations

As with any exercise in experimental evolution used to investigate an ecological phenomenon, there are many limitations that can affect our interpretations and conclusions. The chief concern is that measurements made in the laboratory may not be relevant to organisms in their natural context. However contrived, evolutionary and ecological phenomena still apply to our experimental microcosm, which makes it nevertheless useful for examining ecological questions. Consequently, we are comfortable interpreting our results and discussing them in light of broader ecological theory. Next, our study did not measure changes in phage fitness, but rather, phage productivity. Though this trait can be closely associated with phage fitness and has been previously used as a proxy for fitness (Morley et al., 2015; Turner, Draghi, & Wilpiszeski, 2012), it is clear from our findings that assuming this equivalence is not always appropriate in our system. Importantly, this suggests that the productivity to fitness equivalence might not exist for other virus–host systems. Lastly, our interpretations should consider the potentially complicating role of phenotypic plasticity in these results, as phage phenotypes may be moderated by the host encountered. For example, previous work in T4 has shown that expression of T4 mutant phenotypes can be moderated by bacterial host strain (Benzer, 1957).

Evolutionary genomics data reveals evidence for selection

We observed genetic signatures that are suggestive of selection shaping the evolved populations during niche breadth evolution. Specifically, the high dN/dS ratio, along with the observation that few nonsynonymous mutations reached fixation, indicates selection. The genotypic evidence is suggestive of positive selection as the source of elevated rate of nonsynonymous mutations in structural genes (for a more detailed discussion of other plausible causes, such as mutation bias, see the Appendix A). Though mutation bias can influence the nature of mutations that arrive, newly introduced mutations are still subject to elimination or fixation through the processes of drift and selection. As such, it is probable that fixed (or high frequency) nonsynonymous mutations are of adaptive significance (even if their presence was influenced by mutation bias).

The role of structural genes in bacteriophage niche breadth evolution

Our findings demonstrate that structural genes play a central role in niche breadth evolution: i) The overwhelming majority (82%) of new mutations occurred in structural genes; ii) nearly all (94%) high frequency mutations (>50%) occurred in structural genes; and iii) we detected a significantly higher mutation rate in structural genes. Notably, these rates did not differ across evolutionary histories, which indicates that structural genes are important niche breadth evolution, regardless of the host(s) encountered and level of temporal variability of the environment. Structural genes encode for four different types of virion proteins: capsid, neck, tail, and tail fibers (Miller et al., 2003). Of the 142 new mutations detected in structural genes, 44% occurred in genes coding for tail fiber proteins. We observed mutations in gene 37, which encodes the distal subunit of the LTF (gp37) that directly interacts with bacterial receptors (Bartual et al., 2010). This supports previous evidence suggesting that mutation in gene 37 can alter the niche breadth of T4 (Tétart, Repoila, Monod, & Krisch, 1996). Interestingly, the majority of mutations (74%) detected in LTF genes occurred in the genes encoding the LTF proximal subunits (i.e., gp34 and gp35). This suggests that structures that are only indirectly responsible for host recognition are also important to niche breadth evolution. Mutations in tail genes accounted for 30% of structural mutations; the great majority of those mutations (82%) occurred in genes that code for various baseplate components (e.g., genes 5–10, 27, 48, and 54). Mutations in baseplate genes occurred across all three evolutionary histories, in nearly every population (14 out of 15). This suggests that the baseplate, which regulates infection (Yap et al., 2016), may serve an important role in the evolution of niche breadth. The remaining structural mutations were detected in six genes that code for capsid proteins, with the vast majority (90%) occurring in two of those six genes: soc (small outer capsid) and hoc (head outer capsid). This suggests that soc and hoc proteins, which provide stability to the capsid (Fokine et al., 2004), may also be important to niche breadth evolution. For those interested in a more detailed discussion of mutations, for example, parallel substitutions or distribution across populations, please refer to Appendix A.

CONCLUSIONS

In this study, we used experimental evolution of Escherichia virus T4 to characterize the phenotypic and molecular changes associated with niche breadth evolution. Our findings indicate that populations evolved measurable and meaningful changes in phage productivity. These phenotypic results confirmed some and contradicted other theoretical expectations, which demonstrates the complexity of niche breadth evolution, and the dubiousness of theoretical expectations: Whether generalism truly comes at the expense of high productivity in individual niches is likely context and system‐specific.

Genomic sequencing of evolved populations enabled us to detect a significantly higher rate of mutation in structural genes. Notably, mutation rates of structural genes were similar across evolutionary histories, which indicates that selection in structural genes is important to niche breadth evolution, regardless of conditions that promoted a particular ecological strategy (specialism or generalism). Our study presents compelling evidence that structural genes serve an important role in how phage evolve different niche breadth strategies. This finding further supports the importance of phage structure in the evolution of niche breadth, a result that may have implications for public health surveillance and in therapy. For example, insofar as variation in structural genes is associated with different niche breadth exploitation, perhaps these genes can be manipulated to facilitate a phage's ability to fight particular infections, as in phage therapy (Kortright, Chan, Koff, & Turner, 2019). In addition, if phage genes are associated with expansion of host range, then they may be candidates for a “watch list” approach in disease surveillance, where natural populations are surveyed for mutations associated with exploiting new hosts (Patel, Quates, Johnson, & Ytreberg, 2019).

CONFLICT OF INTERESTS

None declared.

AUTHOR CONTRIBUTIONS

Conceptualization: Jenny Y. Pham, C. Brandon Ogbunugafor, and Daniel L. Hartl. Formal analysis: Jenny Y. Pham and Alex N. Nguyen Ba. Funding acquisition: Jenny Y. Pham, C. Brandon Ogbunugafor, and Daniel L. Hartl. Writing—original draft preparation: Jenny Y. Pham, C. Brandon Ogbunugafor, and Daniel L. Hartl. Writing—review and editing: Jenny Y. Pham, Alex N. Nguyen Ba, C. Brandon Ogbunugafor, and Daniel L. Hartl.

ETHICS STATEMENT

None required.