Temperate Bacteriophages from Chronic Pseudomonas aeruginosa Lung Infections Show Disease-Specific Changes in Host Range and Modulate Antimicrobial Susceptibility

Pseudomonas aeruginosa is a key opportunistic respiratory pathogen in patients with cystic fibrosis and non-cystic fibrosis bronchiectasis. The genomes of these pathogens are enriched with mobile genetic elements including diverse temperate phages. While the temperate phages of the Liverpool epidemic strain have been shown to be active in the human lung and enhance fitness in a rat lung infection model, little is known about their mobilization more broadly across P. aeruginosa in chronic respiratory infection. Using a novel metagenomic approach, we identified eight groups of temperate phages that were mobilized from 94 clinical P. aeruginosa isolates. Temperate phages from P. aeruginosa isolated from more advanced disease showed high infectivity rates across a wide range of P. aeruginosa genotypes. Furthermore, we showed that multiple phages altered the susceptibility of PAO1 to antibiotics at subinhibitory concentrations.

IMPORTANCE Pseudomonas aeruginosa is a key opportunistic respiratory pathogen in patients with cystic fibrosis and non-cystic fibrosis bronchiectasis. The genomes of these pathogens are enriched with mobile genetic elements including diverse temperate phages. While the temperate phages of the Liverpool epidemic strain have been shown to be active in the human lung and enhance fitness in a rat lung infection model, little is known about their mobilization more broadly across P. aeruginosa in chronic respiratory infection. Using a novel metagenomic approach, we identified eight groups of temperate phages that were mobilized from 94 clinical P. aeruginosa isolates. Temperate phages from P. aeruginosa isolated from more advanced disease showed high infectivity rates across a wide range of P. aeruginosa C ystic fibrosis (CF) and bronchiectasis (BR) are respiratory diseases that predispose the individual to chronic infection with opportunistic bacterial pathogens. Pseudomonas aeruginosa is a key pathogen in CF and BR that is associated with increased mortality (1,2). P. aeruginosa can grow as a biofilm, which aids colonization and likely worsens outcomes in the chronically infected lung (3). The difficulty in effective delivery of antibiotics to the lung additionally compromises this treatment, and poor antimicrobial pharmacokinetics may contribute to the development of antimicrobial resistance (AMR) (1,4). P. aeruginosa has a large and flexible genome including a complex accessory genome (5,6). Comparative genomics shows that P. aeruginosa genomes contain a diversity of phage elements that represent a rich source of new genetic material, including antibiotic resistance genes that may be present and undergo frequent recombination (6)(7)(8). Therefore, temperate phages integrated into the bacterial host genome as a lysogen are likely to play an important roles in P. aeruginosa evolution as agents facilitating horizontal gene transfer (HGT) of cargo DNA and additional gene functions. However, changes in integrated phage genomes lead to deactivation of phage induction, and they are therefore no longer capable of lysis or onward HGT (9). Remaining phage gene regions may be a rich source of important gene function and a reservoir for mobilizing phage recombination. The fraction of integrated temperate phages in P. aeruginosa genomes that remain inducible is unknown, since this is rarely tested in comparative genomics studies, limiting our understanding of the potential for phage mobilization in chronic lung infections.
Previous work on the Liverpool epidemic strain (LES) of P. aeruginosa, which was isolated from a CF patient, has shown that this strain harbors multiple inducible temperate phages in its genome (10). These phages remain active and capable of lysis even after many years of chronic CF lung infection, suggesting that they may play an important role in the CF lung (11). Experimental infections in both insect and mammalian hosts demonstrate that the LES temperate phages increase the fitness of P. aeruginosa lysogens in vivo, through phage-mediated killing of competing P. aeruginosa strains (12,13). Furthermore, P. aeruginosa is known to adapt to the lung environment during chronic infection leading to loss of motility (14), mucoid phenotype (15), and biofilm formation (16,17). The LES temperate phages have been shown to contribute to this P. aeruginosa evolutionary adaptation in a sputum-like medium by causing phage insertion-mediated beneficial mutations, particularly in motility and quorum-sensing-associated genes (18). Given these potentially important contributions to P. aeruginosa, we need to understand the potential for phage mobilization across a wider range of clinical P. aeruginosa strains.
Separating individual phage genomes from meta-data is challenging, as related phages share genomic similarity, which increases the probability of assembling chimeric phages. This currently means that it is difficult to characterize phage species at the genome resolution in a mixed population, which requires isolation and propagation of individual phages before genomic comparison and characterization can occur. This can be complicated further by the lack of a receptive bacterial host. We therefore used a novel metagenomic approach and k-mer abundance strategy to resolve individual temperate phages that are mobilizing in complexes from P. aeruginosa genomes. This approach uses the differences in numbers of released phages from the cell during chemical induction and targets these differences in copy number of phage genomes. Here, we characterized the diversity of temperate phages induced from a large panel of P. aeruginosa isolates isolated from patients with either CF or BR. The mobilizable phages were compared at the genomic level, and infectivity profiles were used to determine the similarity and differences between CF and BR. Furthermore, we evalu-ated the effects of temperate phage infection on the antimicrobial susceptibility of P. aeruginosa in vitro.

RESULTS
Phage infectivity profiles. Lambdoid-like phages were chemically induced from each P. aeruginosa isolate using the fluoroquinolone antibiotic norfloxacin. Importantly, all the clinical P. aeruginosa isolates in this panel contained inducible prophages, which is a higher proportion than previously reported in other Gram-negative bacteria such as Escherichia coli (19). Binary networks of phage lysate infectivity against all P. aeruginosa isolates revealed that CF phages are more broadly infectious than BR phages (i.e., higher connectance) ( Fig. 1A and C). However, interestingly, binary networks of BR phage lysates (n ϭ 47) had a more nested structure with higher connectance against mucoid BR isolates (n ϭ 22) ( Fig. 1B and C), suggesting that these phages were better able to infect mucoid phenotypes than nonmucoid P. aeruginosa strains. Semisupervised data reduction modeling (PLS-DA) suggested that the phage lysate infectivity range against P. aeruginosa was greater for phages induced from P. aeruginosa isolated from more advanced stages of disease in both CF and BR (see Fig. S1 in the supplemental material). This allowed stratification of the phage lysates into "early" and "late" disease stages (CF patients were divided into those under the age 16 years or over age 16; BR patients were divided into those less than 10 years since diagnosis or more than 10 years since diagnosis) for subsequent analyses.
Isolating single phage genomes from mixed phage populations induced from a clonal bacterial culture using k-mer count and abundance. We used the biology of temperate phage life cycles and the presumed differences in burst size and diversity to enable 105 complete individual phage genomes to be resolved from the DNA sequencing data. With various numbers of phages in each mixed lysate, a method was implemented to transform the data to distinguish between copy numbers of different viral genomes to resolve individual phages. The Khmer analysis toolkit (20) uses k-mer count and abundance to separate data and is normally used to determine the error k-mer that hinders genome assembly. We modified the python script (calc_median_distribution.py) in Khmer to allow extraction of individual phage sequences that could be assembled into single nonchimeric genomes (see the complete method in Text S1 in the supplemental material) using k-mer distribution. The workflow used to assemble each phage is summarized in Fig. S2.
Analysis of P. aeruginosa phages in CF and BR. Each of the 105 phages was assigned to one of 8 groups (A to H) using a genome comparison dot plot method, Gepard (21) (see Fig. S3). The phage genomes were ordered based on genome similarity and disease etiology. Figure 2 shows whole-genome phylogenetic comparisons aligned with MAFFT, using the L-INS-i algorithm for higher accuracy, which also shows phylogenetic clustering into 8 groups. Basic local alignment (BLAST) was used to determine the closest relative within each phylogenetic clade which was labeled accordingly.
A key aim of this study was to report individual phage genomes; therefore, we only included single contiguous sequence assemblies in downstream analyses (see Table S1). The stringency of our selection criteria removed any phages with fragmented assemblies or low sequence coverage to allow single phage genome comparisons (for partial phage genomes identified within the data, see Table S1). Taxonomic assignment of assembled genomes indicated that Caudovirales phages were enriched in P. aeruginosa isolates in CF and BR (Fig. 2). Groups A to D, F, and G were related to the Siphoviridae family (n ϭ 91), group E was related to the Podoviridae family (n ϭ 13), and group H to the Myoviridae family (n ϭ 1). The distributions of phages between the groups (A to H) and disease etiologies are shown Table S2.
Tables S1 and S2 show both complete and partial genomes determined from the data. We identified phage group A (F10 like) as the most frequently found group (n ϭ 32) in both disease backgrounds. Group E phages (H66/F116 like) were commonly identified in BR. In contrast, group D (D3112-like), C (Phi 297-like), and group G (B3-like) phages were more frequently found in CF. The Phi297-like phages show genome homology to D3 and F116, and this is shown in the Gepard plot (Fig. S3) (22). Interestingly group D phages (D3112 like) (Mu-like transposable phages) (23) and group C (Phi297 like) were found in the same P. aeruginosa metaviromes, suggesting a FIG 1 Shows nestedness and connectance plots demonstrating that temperate phages from P. aeruginosa isolated from patients with advanced respiratory infection can infect a wider range of P. aeruginosa genotypes. The nestedness and connectance plots map the cross-infection data from mixed temperate phage communities induced from clinical P. aeruginosa isolates. (A) The binary nested network for the complete 94 phage lysates against 94 P. aeruginosa isolates, ordered by nestedness within equally sized quadrants (47 CF phages versus 47 CF P. aeruginosa isolates, 47 BR phages versus 47 CF P. aeruginosa isolates, 47 CF phages versus 47 BR P. aeruginosa isolates, and 47 BR phages versus 47 BR P. aeruginosa isolates). Each black square represents an individual interaction (infection event) between one phage strain (x axis) and one bacterial strain (y axis). White squares represent lack of infection. Note the CF-derived phages are capable of infecting BR P. aeruginosa isolates more frequently than BR phages can infect CF P. aeruginosa, illustrated through higher values for connectance (0.693). (B) The binary nested network representation of infection of the mucoid BR P. aeruginosa isolates (n ϭ 22) with the entire cohort of mixed temperate phage communities induced from BR P. aeruginosa isolates (n ϭ 47). (C) Nestedness and connectance values for the results shown in panels A and B. possible co-lineage. We also identified a cytotoxin-converting PhiCTX-like phage (group H), which carries the CTX gene; this phage was found in a single BR P. aeruginosa isolate (24).
Antimicrobial susceptibility of twenty PAO1 lysogens. From each clinical group (CF, Ͻ16 years of age, Ͼ16 years; BR, Ͻ10 and Ͼ10 years since diagnosis), 5 phage lysates were randomly selected to create P. aeruginosa PAO1 lysogens. Of the 20 lysogens, 8 showed altered antimicrobial susceptibility to clinically relevant antimicrobials compared to that of wild-type PAO1 when testing a concentration subinhibitory to the MIC for PAO1. Figure 3 shows the altered phage-mediated growth of PAO1 in the presence/absence of ceftazidime (Fig. 3A), colistin (Fig. 3B), meropenem (Fig. 3C), and piperacillin ( Fig. 3D) (25). Seven lysogens significantly decreased the inhibition of growth to meropenem and piperacillin; one significantly increased inhibition of growth FIG 2 Legend (Continued) evolutionary history was inferred using the neighbor-joining method, with the sum of branch length of the optimal tree of 3.04826762. The percentages of replicate trees in which the associated taxa clustered together in the bootstrap test (1,000 replicates) are shown next to the branches; the tree was rooted at the midpoint. The evolutionary distances were computed using the p-distance method and are in the units of the number of base differences per site. All ambiguous positions were removed for each sequence pair. There were a total of 129,881 positions in the data set. Evolutionary analyses were conducted in MEGA7.

FIG 3
Antimicrobial susceptibility of lysogens to clinically relevant antibiotics. The data shown are after 9 h of incubation: ceftazidime, 0.08 g/ml (A); colistin, 1.6 g/ml (B); meropenem, 0.08 g/ml (C); and piperacillin, 0.8 g/ml (D). *, P Յ 0.05 by one-way ANOVA with Dunnet's post hoc test for parametric variables or a one-way ANOVA Kruskal-Wallis test with Dunn's multiple-comparison test. The lysogens were grouped (and color coded) based on disease etiologies. to meropenem. Antibiotic sensitivity was also compared to disease etiology, with only childhood CF phage PAO1 lysogens showing significant reduction in sensitivity to meropenem and piperacillin (see Fig. S4). The genomes of the 20 lysogens were sequenced and compared to that of wild-type PAO1 using Mauve (26) to determine phage integration sites (see Table S3). Importantly, we found that none of the phages integrated within putative genes that could be associated with this phenotype. The specific phage integrated into each site in the PAO1 chromosome is also detailed in Table S3.
Phage genome search for AMR gene prediction. The 105 phage genomes were searched using 2 databases (CARD [27] and ARDB [28]) for known and previously reported AMR genes; none were identified in the phage genomes.

DISCUSSION
This is the largest pan-genome study of P. aeruginosa temperate phages, which mapped the interactions between phage and bacteria and their broad infectivity across chronic respiratory disease-related P. aeruginosa isolates. We investigated a broad range of P. aeruginosa isolates from patients at different stages of disease, and to our knowledge, this is the largest reported cross-infection study of lysogenic phages undertaken for clinically derived P. aeruginosa isolates, with over 8,800 infections. We first showed that P. aeruginosa isolates from these diverse clinical sources have high levels of inducible phages compared to those of other Gram-negative bacteria (19). This illustrates a significant reservoir of mobile DNA that could play a major role in bacterial adaptation and evolution in the lung. This finding combined with a recent study by Nguyen and colleagues (29) demonstrating that phages are able to interact and translocate across epithelial surfaces adds a further layer of largely unexplored complexity to chronic respiratory infections.
We found a positive correlation between phage infectivity and disease progression in both CF and BR, which supports a study by James and colleagues that followed LES phage populations in LES-infected CF patients over a period of 28 months (11). Our approach here is more representative of phage induction in the lung, as it utilizes the infectivity of all the temperate phages, mobilized from each P. aeruginosa genome under selective pressure. This provides a panoramic snapshot of phage mobilization in vivo that may illustrate the broad reach of HGT in chronic respiratory disease.
One of our key findings is the hierarchical infectivity of CF-and BR-derived phage lysates by bipartite ecological network modeling. We saw that chronicity of infection was associated with greater levels of infectivity. These infectivity profiles supervised by age and time in CF and BR, respectively, show significant differences that allowed us to stratify this cross-sectional panel into 2 groups per disease. Previous studies have demonstrated that P. aeruginosa isolated from CF patients can cross-infect patients already infected with other P. aeruginosa (30)(31)(32), in particular, the Liverpool epidemic strain (LES). We have previously shown using Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of these phage metagenomes that the size of the phage accessory genome increased with the progression of disease. Furthermore, we identified gene markers that were associated with survival in CF or BR patients with chronic P. aeruginosa infections (33). Through classifying these phage-host interactions, we showed the breadth of HGT that would occur on phage induction where we have previously described the function that can be disseminated directly to the bacterial host from the phage metagenomes (33). Similar studies have found that temperate phages mediate pathogen fitness in the lungs of CF patients and demonstrate their clinical importance in the pathophysiology of the disease (10,11,18). The broad infectivity of phage lysates derived from adult CF patients when infecting P. aeruginosa from other clinical sources suggests there are additional risks beyond our current understanding of "bacterial cross-infection," where phage-mediated HGT adds another level of complexity.
Importantly, we developed a novel method for resolving individual phage genomes from metaviromes by modifying a program (Khmer [20]) used to remove the error k-mer from genome sequencing data. This method has allowed resolution of a variety of phage genotypes present in chronic respiratory disease-related P. aeruginosa. Given the importance of transposable phages in driving adaptive evolution of bacterial pathogens within hosts, the findings of this study highlight a number of temperate phages to focus on in prospective studies. For example, the group A or F10-like phages were the most commonly found and more diverse group (n ϭ 32) across this panel of P. aeruginosa; however, despite their prevalence, there is currently limited research on the biology of these temperate phage (34). Other noteworthy examples include F116like phages (group E); these have been shown to easily exchange bacterial chromosome genes (35), suggesting it may have a role in bacterial adaptation in the chronically infected lung. Phi297 has been described as producing halo-morphologies and has a narrow host range (22); in this study, the Phi297-like phages were predominantly found in CF isolates and more likely to be associated with a broader infectivity profile.
The phage genomes from group D (D3112 like) and group C (Phi297 like) in our cohort were intriguingly detected together, a marker that in tandem, they are potentially ecologically important, displaying a symbiotic relationship. We tested this further, and even when full contiguous sequences of one phages were present, the secondary phage was always found within the other k-mer peak assemblies.
Kim et al. showed that if P. aeruginosa develops resistance to PA1Ø phage, which shares genome homology to D3112 (a pilus-infecting phage), then biofilm formation was reduced, which leads to increased antibiotic susceptibility (36). This provides evidence that suggests that cocarriage of phages is required and may constitute an evolutionary advantage, something that has not been previously shown in P. aeruginosa. It is noteworthy that even though each clade relates to a previously described phage infecting P. aeruginosa, the genomic diversity within each phage group is extremely high and requires further study. Current research focuses on identifying genes in bacterial backgrounds and studying the incidence of AMR, which is a key problem in P. aeruginosa infection in CF and BR. In this study, we demonstrate that phages modulate the growth of the lab strain PAO1 in the presence of antimicrobials and in the absence of known resistance genes. We show that 7 of the 20 phages decreased the susceptibility of P. aeruginosa lab strain PAO1 against one or more of the clinically relevant antipseudomonals and that 1 increased the susceptibility of P. aeruginosa at concentrations below the clinical MIC for PAO1; yet, the concentration of the antibiotic inhibited the growth of PAO1. Previous work by Fothergill et al. illustrated differential induction rates and, in some cases, reduction of phage release when LES P. aeruginosa isolates were challenged with commonly prescribed antibiotics (37). Microbiome studies of the gut have shown that antibiotic treatment enriched the phage metagenome and increased transfer of AMR genes under conditions of physiological stress (38,39). Importantly we found no gene(s) associated with these small but significant changes in bacterial susceptibility to antibiotics using curated CARD and ARDB database searching. We therefore hypothesize that this association is mediated by phages and that they subvert normal cellular and metabolic pathways of P. aeruginosa. This was previously described for phages transferring genes horizontally in Escherichia coli (40). Here, we showed that phages 24 (CF) and 200 (BR) when integrated into PAO1 as lysogens have the greatest impact reducing the susceptibility to the antibiotics tested. We propose that a decrease in susceptibility may increase tolerance, which in turn may progress to AMR (41). Phage 24 is group A (F10 like) and the most heterogeneous compared to the other F10-like phages isolated here. Similarly, phage 200 is group D (D3112) and integrates at a different point in the genome than other group D phages. Observationally, phages associated with earlier stages of disease in CF provide P. aeruginosa with an increased resistance to antibiotic therapy. Does this mean that these phages are supportive in early colonization of the CF lung?
To date, there are no treatment algorithms sufficient to meet the demands of these new clinical challenges. The ability of phages to subvert the physiology of P. aeruginosa is worrying, especially given the lack of an associated gene target available for sequence output file was directly pipelined into the Khmer toolkit v1.1, which was used for its designed ability to remove very low-level bacterial contamination and erroneous k-mers from the viral sequence data. The Khmer count table and k-mer abundance histogram clusters low abundance and poor sequence data that would be linked to any residual bacterial chromosomal DNA (20). The Khmer generated k-mer abundance data were graphed in Excel, and the error k-mer peak was manually selected and removed. The calc-median-distribution.py script part of Khmer (20) was modified and renamed k-mer_extraction.py (https://github.com/rmadnantariq/k-mer_extraction). The script was altered to select out the k-mers associated with an abundance peak distribution described in detail in the methods in the supplemental material.
Nucleotide dot plot analysis using Gepard The 105 phage genomes were concatenated together based on genome similarity and ordered based on disease etiologies. The dot plot was generated with Gepard (21) using a word size of 10 and sliding window of 25.
Lysogen assembly. CLC genomic workbench v11 was used to generate the assembly of the naive PAO1 lab strain. This assembly was used as a reference in Ragout post de novo assembly of the lysogens using SPAdes (v3.5.0). The assemblies were annotated using Prokka v1.22. The naive PAO1 and referenceassembled PAO1 lysogens were compared using Mauve to track the phage and their insertion sites.
Data availability. All data generated or analyzed during this study are included in the published article (and its supplemental material files). The raw fastq files were deposited in GenBank with the BioProject accession PRJNA503342. Resolved 105 phage genome accession numbers can be found in Table S1.

ACKNOWLEDGMENTS
We thank the research nurses and clinical staff at the Freeman Hospital and Royal Victoria Infirmary, Newcastle upon Tyne, UK, for all their help collecting the samples that yielded the P. aeruginosa isolates used in this investigation. We also thank all of the patients that kindly participated in this study. Our thanks to the Titus Brown and Khmer developers group for allowing us to edit the bioinformatics script. Finally, we acknowledge Nathan Brown for his introduction to PriceTI and providing bioinformatic support.