Selective Sweeps and Parallel Pathoadaptation Drive Pseudomonas aeruginosa Evolution in the Cystic Fibrosis Lung.

UNLABELLED: Pulmonary infections caused by Pseudomonas aeruginosa are a recalcitrant problem in cystic fibrosis (CF) patients. While the clinical implications and long-term evolutionary patterns of these infections are well studied, we know little about the short-term population dynamics that enable this pathogen to persist despite aggressive antimicrobial therapy. Here, we describe a short-term population genomic analysis of 233 P. aeruginosa isolates collected from 12 sputum specimens obtained over a 1-year period from a single patient. Whole-genome sequencing and antimicrobial susceptibility profiling identified the expansion of two clonal lineages. The first lineage originated from the coalescence of the entire sample less than 3 years before the end of the study and gave rise to a high-diversity ancestral population. The second expansion occurred 2 years later and gave rise to a derived population with a strong signal of positive selection. These events show characteristics consistent with recurrent selective sweeps. While we cannot identify the specific mutations responsible for the origins of the clonal lineages, we find that the majority of mutations occur in loci previously associated with virulence and resistance. Additionally, approximately one-third of all mutations occur in loci that are mutated multiple times, highlighting the importance of parallel pathoadaptation. One such locus is the gene encoding penicillin-binding protein 3, which received three independent mutations. Our functional analysis of these alleles shows that they provide differential fitness benefits dependent on the antibiotic under selection. These data reveal that bacterial populations can undergo extensive and dramatic changes that are not revealed by lower-resolution analyses. IMPORTANCE: Pseudomonas aeruginosa is a bacterial opportunistic pathogen responsible for significant morbidity and mortality in cystic fibrosis (CF) patients. Once it has colonized the lung in CF, it is highly resilient and rarely eradicated. This study presents a deep sampling examination of the fine-scale evolutionary dynamics of P. aeruginosa in the lungs of a chronically infected CF patient. We show that diversity of P. aeruginosa is driven by recurrent clonal emergence and expansion within this patient and identify potential adaptive variants associated with these events. This high-resolution sequencing strategy thus reveals important intraspecies dynamics that explain a clinically important phenomenon not evident at a lower-resolution analysis of community structure.

INTRODUCTION

Cystic fibrosis (CF) is the most common fatal genetic disease among Caucasians (1, 2) and is characterized by chronic polymicrobial lung infections that are frequently dominated by opportunistic pathogens, such as Pseudomonas aeruginosa and species from the Burkholderia cepacia complex (3, 4). Polymicrobial communities in adult CF airways appear to be highly resilient and resistant to change, despite the administration of sequential antibiotic regimens targeting these dominant pathogens (5, 6). However, our perception that these communities are largely static is predominantly driven by culture-dependent studies of sputum and bronchoalveolar lavage fluid, microbiome studies using subregions of the 16S rRNA gene, and long-term longitudinal comparative genomics studies (7–9). While these approaches have tremendously informed our understanding of CF lung microbiology, they have not been performed at the temporal and taxonomic resolution required to characterize short-term intraspecies population dynamics. Since selection on microbial communities acts at the level of individual clones over the short temporal scale of antimicrobial treatments and lung disease exacerbations (10), a higher-resolution taxonomic and temporal approach is required to identify, understand, and ultimately manipulate the factors that control the composition and function of polymicrobial communities.

CF lung microbial communities persist in the face of powerful and complex pressures imposed by episodic exacerbations of lung inflammation and courses of antimicrobials, which typically include both long-term treatment with maintenance antibiotics as well as short-term, recurrent courses of other antimicrobials administered alone or in combination (11, 12). A number of studies have used comparative genomics to address how CF pathogens and P. aeruginosa in particular adapt to this dynamic lung environment. These have shown that airway-adapted isolates display very high levels of phenotypic diversity, strong mutational parallelism (recurrent mutations in the same genes among independent lineages), the emergence of hypermutable strains, mutations in global regulatory genes, and long-term persistence of a small number of clonal lineages (13–29). While these studies have furthered our understanding of the long-term evolutionary dynamics of P. aeruginosa in the CF lung, the large majority have been based on isolates collected several years apart—in some cases, 35 or more years apart—and none have focused on the immediate impact of antimicrobial therapies on the overall population of P. aeruginosa in the CF lung over a very short time span (19, 21, 24, 26, 27). This is a critical gap in our knowledge since clinical assessments and antimicrobial therapies are often modified as frequently as weekly (30).

A fundamental assumption of antimicrobial therapies is that they select against and thereby drive down the population size of target pathogens. A nearly unavoidable consequence of these therapies is the selection for resistant clones of the pathogen, which may ultimately dominate the population. The application of a selective pressure to microbial populations commonly leads to selective sweeps, which involve the spread of a clone carrying an adaptive mutation through the population and subsequent displacement of other preexisting and now less fit clones (31–34). Selective sweeps will periodically reoccur in a population as new beneficial mutations arise and new selected clones replace ancestral clones. This pattern of recurrent selective sweeps has historically been called periodic selection (35). While selective sweeps have been well documented in microbial populations (particularly with respect to viral populations, as reviewed in reference 36), the presence and significance of periodic selection in human health and infectious diseases are not known (37).

Clonal evolutionary dynamics such as selective sweeps have been largely overlooked in microbiome studies due to the lack of resolution and power. For example, 16S rRNA gene sequencing can only reliably resolve taxa to the genus level (38), while full metagenome sequencing suffers from the inability to phase genetic variation (assign variants into genotypes) (39). These limitations have restricted our ability to observe clonal evolutionary dynamics, identify their drivers, and assess their resulting impact on the overall microbial community as well as the host.

In the case of CF-associated lung infections, longitudinal population sampling (e.g., multiple clones assessed from each clinical specimen) performed over the short term will reveal the impact of intensive antimicrobial therapies on the dynamics of the target pathogen populations. Understanding intraspecies pathogen evolution over this time scale is essential to assess whether current treatment strategies effect desirable change of pathogen populations or if alternative treatment strategies, coupled with a better understanding of pathogen evolutionary potential, could lead to improved treatment outcomes.

In this study, we sequenced the genomes of 233 P. aeruginosa isolates from a single CF patient collected from 12 specimens over a 1-year period to assess very-short-term evolutionary dynamics. We identify the occurrence of a very recent common ancestor for the overall population, the emergence of multiple clonal lineages within a population that has remained relatively constant in size, and strong signals of positive selection—genomic features consistent with recurrent selective sweeps or periodic selection. We also find a very strong signal of parallel pathoadaptation at loci associated with virulence and antibiotic resistance and show how some of these variants confer fitness benefits that are dependent on the environment.

RESULTS

This study focuses on patient CF67, a 34-year-old female CF patient (of CF transmembrane conductance regulator [CFTR] genotype ΔF508/del2–3) who had been chronically infected with P. aeruginosa for at least 12 years prior to the initial study sample. The patient had advanced lung disease with a forced expiratory volume in 1 s (FEV1) of 21% of the predicted value throughout the course of the study. She provided 12 expectorated sputum specimens for evaluation over an ~1-year period (with the first specimen collected 350 days prior to the last specimen) (Fig. 1). The specimens were not collected uniformly over the course of the study, with 11 being collected within the last 126 days of the study.

FIG 1 

Antibiotic treatment history and the relative abundance of the two clades over time. (A) Black bars indicate antibiotic administration, and hashed bars indicate intermittent exposure in that time block. The method of antibiotic administration is shown as intravenous (iv), inhaled (inh), or oral (po). Sputum samples were collected at the time points indicated by green lines, which extend to panels B and C. (B) Relative abundance at the genus level is shown as a percentage. (C) Relative abundance of clade A (blue shades) and clade B (red shades) genotypes over time. Genotypes are defined as whole-genome sequences differing by one or more SNPs segregating in at least two isolates. The same color and shading are used to identify specific genotypes across specimens.

During the study period, the patient was prescribed 28 distinct courses of antibiotics, which were administered by three different routes (oral, intravenous, and inhaled) in addition to continuous oral azithromycin (Fig. 1A; see Table S1 in the supplemental material). Despite this complex antimicrobial regimen, microbiome analysis via Illumina sequencing of the V5 to V7 hypervariable regions of the 16S rRNA gene (38) reveals that Pseudomonas remained the dominant genus throughout (Fig. 1B). Quantitative PCR of the 16S rRNA gene was used to quantify the bacterial load from each of the twelve specimens, revealing no appreciable change in the absolute level of bacteria recovered from each specimen over time (mean ± standard deviation [SD], 0.395 ± 0.185 ng/ml; regression slope, −0.0005x).

We collected multiple P. aeruginosa isolates from each of the 12 specimens (range, 18 to 20 per specimen) in order to obtain a finer-resolution analysis of the population-level dynamics. An effort was made to maintain overall morphotype relative abundance in our sampling as previously described (40).

Genomic diversity and phylogenetic analyses identify two genetically distinct clades.

A total of 233 P. aeruginosa isolates underwent whole-genome sequencing (7 of the total of 240 isolates were removed from the analysis due to quality control issues), yielding a median coverage depth of 54× (see Fig. S1 in the supplemental material). The genome of one isolate was assembled de novo for use as a reference for assembly of the remaining isolates. This assembly had an N50 of 329,133 bp and 46 contigs, with the smallest contig being 1,170 bp and a total genome size of 6,367,116 bp. We used a conservative pipeline to call single nucleotide polymorphisms (SNPs) that required the agreement of three algorithmically distinct reference assembly methods and then reevaluated each SNP position using a more relaxed inspection of the pileup data (see Fig. S2 in the supplemental material). The pipeline identified 107 segregating SNPs among the 233 isolates, of which 54 were phylogenetically informative (SNPs segregating in at least two isolates), while 39 were found in isolates collected from multiple time points (multispecimen SNPs). Additionally, we identified insertion/deletion (indel) variants by realigning the sequencing reads to initial candidates and then reviewing all putative indel positions in all isolates. This process revealed 58 indels, of which 28 were found in at least two isolates from different collection time points. Sequencing of the entire population of all P. aeruginosa isolates bulk extracted from one specimen did not reveal any mutations that were not present in the individual isolate sequencing, which indicated that our sampling of diversity is close to saturation (see Fig. S3 in the supplemental material).

The data were analyzed by both network-based (neighbor-net) and traditional (maximum likelihood) phylogenetic methods (Fig. 2). For the latter, the position of the root is very similar, regardless of whether it is identified by incorporating P. aeruginosa LESB58 into the analysis, selecting it based on the position of the strains isolated in specimen 1, or through midpoint rooting. Visual inspection of both trees revealed an unusual topology. Approximately one-third of the isolates share more structured relationships with relatively long branches, and these are distinctly separated from another set of isolates that have relatively short branches originating from a very recent common ancestor. Based on this observation, we designate two clades: clades A and B. While the clades are discrete and distinct in the network-based analyses, the rooted phylogenetic analysis shows that clade B is nested inside clade A, making clade A technically paraphyletic since it is not monophyletic without the clade B strains. Despite this, we will refer to the two groups as clades for simplicity. Importantly, all analyses comparing the two clades are performed on discrete sets of strains, where all strains are designated as belonging to either clade A or clade B as determined by the network-based phylogenetic analysis.

FIG 2

Maximum likelihood and network-based (neighbor-net) phylogenetic analyses. The phylogenetic structure of 233 P. aeruginosa isolates was characterized based on genome-wide single nucleotide polymorphisms (SNPs) using the maximum likelihood (A) and neighbor-net (B) algorithms. The structure of the resulting tree revealed two populations, clades A (blue) and B (red). The 152 isolates in clade B produce a star phylogeny consistent with a recent expansion of a clonal population, while the 81 isolates in clade A show longer branches with more phylogenetic structure. (A) The scale bar shows genetic distance using the maximum composite distance and all segregating SNPs. The strain names shown on the tree are a composite of the specimen number from which the clone was isolated, an arbitrary clone letter, and the clade designation. (B) Individual strain names at the tips of each branch have been replaced with pie charts indicating the distribution of dates during which the strains were sampled (indicated by the circular legend). The scale bar indicates genetic distances.

Clade A has a lower overall relative abundance than clade B (81 isolates versus 152 isolates for clades A and B, respectively [Fig. 1 and Table 1]) but is significantly more genetically diverse, with an average of 4.80 ± 3.97 (mean ± SD) pairwise differences among clade A isolates versus 2.84 ± 1.60 differences for clade B isolates (P <<< 0.0001, t test) (Fig. 3). Neither clade shows a significant recombination signal, with a four-gamete recombination test identifying a minimum of only one recombination event in the entire sample (41). This is supported by the nearly complete lack of reticulation in the neighbor-net phylogenetic analysis. The relatively uniform distance of clade B strains from their common ancestor seen in the phylogenetic analyses is reminiscent of a “star phylogeny” pattern of divergence, which is often associated with recent purges of genetic variation, such as those seen during population expansions or selective sweeps (42). There are six segregating mutations (four SNPs and two indels) that have variants fixed between the two clades. This is significant since fixed differences reflect population subdivision.

TABLE 1 

Genetic diversity

Variable	No. (%)
	Clade A	Clade B	Total
Isolates	81 (34.8)	152 (65.2)	233
SNPs	72 (67.3)	35 (32.7)	107
Indels	31 (53.4)	28 (48.3)	58a
Phylogenetically informative SNPs	36 (64.3)	20 (35.7)	56
Multi-time SNPs	21 (53.8)	18 (46.2)	39 (57.1)
Multi-time indels	13 (46.4)	16 (57.1)	28a
Genetic distance, b avg ± SD	0.052 ± 0.04	0.029 ± 0.03	0.047 ± 0.04
No. of SNP differences, avg ± SD	4.80 ± 3.97	2.84 ± 1.60	6.49 ± 4.41
No. of indel differences, avg ± SD	4.69 ± 3.11	3.05 ± 1.91	5.03 ± 3.09

One indel segregates in clade A and clade B.

Maximum composite likelihood distance.

FIG 3 

Frequency spectrum of SNPs in clades A and B. (A) The frequencies of the number of SNP differences between isolates within clade A and B reveal different distributions. While in clade A, most pairwise SNP differences occurred evenly, in clade A, pairwise SNP differences lower than 5 are overrepresented. (B) Distribution of pairwise SNP differences within clades A and B shows different profiles (P <<< 0.0001, t test).

Time to common ancestry.

We used the Bayesian approach implemented in BEAST (43, 44) to infer the time to most recent common ancestry (tMRCA) for the total population and each clade individually. The total population of 233 isolates coalesces to an MRCA of 993.24 ± 9.09 days (2.72 years) prior to the last specimen, while clade A isolates coalesce to a tMRCA of 956.75 ± 8.93 days (2.70 years) prior to the last specimen, and clade B isolates coalesce to an tMRCA of 329.42 ± 1.80 days (0.90 years) prior to the last specimen. We compared our isolates to other fully sequenced P. aeruginosa isolates (see Fig. S4 in the supplemental material) and found the Liverpool epidemic strain LESB58 to be the closest reference strain, with a tMRCA of 24.06 years prior to the last specimen (45, 46). The probability that our sample of 233 isolates captures the true tMRCA of the population can be calculated by (n − 1)/(n + 1) = 0.991, where n is the sample size (47).

Our analysis supports a population model in which clade A is a longer-lived ancestral population that has accumulated more genetic variation than clade B. Clade B, on the other hand, appears to be a more recently emerged population that arose between collection of the first and second specimens. The relatively recent coalescence of the total sample and clade B in particular suggests that we may have captured two periods of clonal emergence, expansion, and replacement. This is supported by the fact that clade A is more abundant in the earlier specimens, while clade B dominates in the vast majority of the later specimens (Fig. 1).

Selection analysis.

The emergence and spread of a new clone are called a selective sweep when driven by positive selection; therefore, we sought to identify evidence of selective pressure driving the emergence of clade B. Analysis of the 107 segregating sites in the total population gave a ratio of nonsynonymous to synonymous mutations (dN/dS ratio) of 1.01, which is consistent with the null model of neutral evolution but which does not rule out more complex selective patterns (Fig. 4 and Table 2). Importantly, this value is also consistent with a nonequilibrium population, in which the short time frame of the study allows us to capture mutations as they arise but before their fate has been determined by natural selection or genetic drift. Consequently, it is probable that some of the observed SNPs are deleterious to some degree and will be culled from the population in the longer term. To minimize this confounding effect, we elected to focus on only those 39 SNPs that segregated over at least two sampling time points (multispecimen SNPs). It would be expected that removing transient deleterious mutations from the analysis should decrease the relative number of nonsynonymous mutations, thus decreasing dN/dS ratio. Despite this, we find a dN/dS ratio of 1.43, consistent with weak positive selection acting on mutations maintained over multiple time points. When we split the sample based on the inferred clade structure, we find that the clade A, multispecimen segregating SNPs gave a dN/dS ratio of 1.04, compared to a dN/dS ratio of 2.42 for clade B. Thus, clade B is under greater positive selective pressure than clade A. When we further break down the analysis into the 12 individual specimens, we find a very consistent pattern, with clade B maintaining a higher dN/dS ratio relative to clade A throughout, with the exception of specimen 9 (Fig. 4). It is also intriguing the dN/dS ratio observed in clade A starts at 1.3 and then drops to 1.0 or below after specimen 1, once again with the exception of specimen 9. While there are many potential factors that can influence dN/dS ratios, and their calculation is highly dependent on sampling variance, the overall pattern observed in these data supports stronger and more sustained positive selection in clade B than clade A. Notably, the reversal of the dN/dS ratio trend noted in specimen 9 coincided with the cessation of aztreonam during that sampling period.

FIG 4 

dN/dS ratio calculation over time. Shown is the nonsynonymous over synonymous substitution rate (dN/dS) for multispecimen SNPs calculated for each sample time and partitioned based on clade membership. Clade B arose between sampling specimens 1 and 2; therefore, there are no data for clade B.

TABLE 2 

Selection analysis results

Variable	No. (%)			dN/dS ratio
	Nonsynonymous	Synonymous	Intergenic
All SNPs	73 (68.2)	25 (23.4)	9 (8.4)	1.01
Multitime SNPs
All	29 (74.4)	7 (17.9)	3 (7.7)	1.43
Clade A	15 (71.4)	5 (23.8)	1 (4.8)	1.04
Clade B	14 (77.8)	2 (11.1)	2 (11.1)	2.42

The dN/dS results supporting positive selection are further bolstered by the significantly negative result for both a Tajima’s D test (48) and Fu and Li’s test (49). Tajima’s D = to −2.08 (P = 0.002), while Fu and Li’s D* = −5.58 (P = 0.001) and D = to −5.60 (P < 0.001, using the P. aeruginosa LESB58 strain as an outgroup). No significant differences were observed between the two clades in these analyses. These results indicate that there has been a recent population expansion after a selective sweep or population bottleneck. Given that there was no change in the absolute size of the Pseudomonas population in this patient, the best explanation for this expansion is the replacement of one clone by another.

Parallel pathoadaptation.

A striking pattern in the data is the number of loci that were independently mutated multiple times. The occurrence of multiple independent mutations within the same locus can be a very strong indication of adaptive parallel evolution and strong selection on the locus (50–52). Nineteen loci were mutated multiple times (Table 3). Thirty-three of the 107 SNPs (30.8%) occur in these loci, with 15 (38.5%) being multispecimen SNPs, while 20 of the 58 indels (34.5%) are found in loci with multiple independent mutations, 11 of which are multispecimen indels (39.3%).

TABLE 3 

Multiply mutated loci

Encoded protein	Locus	No. of SNPs/indels	Probability a
Murein tripeptide ligase (Mpl)	PA4020	1/5	1.20 × 10−19
Anti-sigma factor (MucA)/sigma factor (AlgU)b	PA0763/PA0762	2/3	1.17 × 10−14
Alginate-C5-mannuronan epimerase (AlgG)	PA3545	5/0	7.31 × 10−16
Conserved hypothetical protein	PA4701	2/2	3.56 × 10−11
Transcriptional regulator (AmpR)/β-lactamase precursor (AmpC)	PA4109/PA4110	4/0	3.56 × 10−11
Penicillin-binding protein 3 (PBP3)	PA4418	3/0	2.70 × 10−08
Multidrug efflux membrane fusion protein (MexA)/multidrug efflux transporter (MexB)	PA0425/PA0426	0/3	1.08 × 10−07
Transcriptional regulator (MexT)	PA2492	2/0	1.64 × 10−04
Outer membrane lipoprotein precursor (OprI)c	PA2853	2/0	1.64 × 10−04
Hypothetical protein	PA3093	2/0	1.64 × 10−04
Hypothetical proteind	NCGM23404	2/0	1.64 × 10−04
Aminopeptidase P (PepP)/ubiquinone biosynthesis protein (UbiH)c	PA5224/PA5223	0/2	3.29 × 10−04
ld-Carboxypeptidase (LdcA)	PA5198	1/1	1.64 × 10−04
Conserved hypothetical protein	PA5133	1/1	1.64 × 10−04
Conserved hypothetical protein/hypothetical protein	PA4962/PA4961	0/2	3.29 × 10−04
Transcriptional regulator (LasR)	PA1430	1/1	1.64 × 10−04
Multidrug resistance operon repressor (MexR)	PA0424	1/1	1.64 × 10−04
Probable transcriptional regulator	PA0535	1/1	1.64 × 10−04
d-Amino acid dehydrogenase, small subunit (DadA)	PA5304	0/2	1.64 × 10−04

Calculated based on the probability of resampling with replacement any locus n times, given a genome size of N. P = (1/N)^(n − 1). We used (n − 1) since we are calculating the probability for any locus, rather than a specific locus.

Mutations that occurred in two adjacent loci and which may be part of a common regulon according to the Pseudomonas Genome Database (97). The probability calculation was adjusted by dividing the total number of loci in the genome in half.

One mutation occurred in the intergenic region flanking this locus.

This locus is not found in PAO1. The homolog with highest similarity is in P. aeruginosa NCGM2.S1 (100).

We compared the mutational class spectrum (relative frequency of nonsynonymous, synonymous, and intergenic changes) in our entire data set to that seen in loci with multiple mutations to determine if the parallel evolution is due to selection or simply a higher local mutation rate. We would expect an elevated mutation rate to maintain the same ratio of nonsynonymous to synonymous changes, while adaptive parallel evolution should predominantly be driven by nonsynonymous changes. We observed a mutational class spectrum of 68.2% nonsynonymous, 23.4% synonymous, and 8.4%, intergenic mutations in our full data set and a spectrum of 91.7%, 0.0%, and 8.3%, respectively, for the loci with multiple mutations. This highly significant skew toward more nonsynonymous mutations (chi-square test, P = 1.5 × 10−7) strongly supports selection-driven parallel evolution.

Many of the loci with multiple mutations are known to play important roles in P. aeruginosa virulence in the setting of CF. Among them is the gene that encodes a murein peptide ligase (Mpl), which has been shown to increase expression of AmpC and resistance to the β-lactam antibiotics piperacillin and ceftazidime (53). This locus carries five independent indels and one SNP. Given a genome size of 6,082 genes and making some simplifying assumptions concerning uniformity of gene size, the probability of any gene being mutated six times is 1.2 × 10−19. Similarly, genes involved in alginate biosynthesis, such as the MucA and AlgU cluster and the locus that encodes AlgG, were each mutated five times (P = 1.17 × 10−14 and P = 7.3 × 10−16, respectively). Many other loci carrying multiple independent mutations also influence antimicrobial resistance, such as the operon encoding multidrug resistance-related loci MexA and MexB, which was mutated four times (P = 1.08 × 10−07), the transcriptional regulator AmpR and its neighboring β-lactamase precursor, AmpC, which were mutated four times (P = 3.6 × 10−11), and the gene encoding penicillin-binding protein 3 (PBP3), which was mutated three times (P = 2.7 × 10−8). Six other loci were mutated twice (P = 1.6 × 10−4 for any one locus to receive two independent mutations), including the gene encoding the transcriptional regulator MexT.

Identification of potential adaptive variants.

To detect possible adaptive variants, we determined the ontology of multispecimen segregating variants. Many loci with segregating variation in our sample are implicated in antibiotic resistance and adaptation to the host (see Table S2 in the supplemental material). The six fixed differences distinguishing the two clades include SNPs in genes encoding the two-component regulator NarX, a phosphomannose isomerase, an arginine/ornithine succinyltransferase, and the fimbrial subunit CupC1. All but the last one are nonsynonymous substitutions. One of the fixed indels leads to a frameshift mutation in the locus that encodes AlgF, which controls the addition of acetyl groups to alginate. The second fixed indel is located in the intergenic region upstream from the gene encoding DNA polymerase I (polA). Beyond these fixed mutations, there are SNPs in the gene encoding a probable nonribosomal peptide synthetase and alginate-C5-mannuronan epimerase AlgG that are also fixed between the two clades, with the exception of a single clade A specimen 1 isolate that carries the clade B variant at both sites. Additionally, there are SNPs in genes encoding a β-lactamase precursor and the transcriptional regulator MexT that are found only in clade B isolates and that first appear in specimen 2, which is the first specimen containing clade B isolates. These SNPs are found in 29.0% and 17.8% of clade B strains, respectively. Finally, an indel located in the gene that encodes a murein tripeptide ligase (Mpl) segregates exclusively in clade B, where it is found in 76.3% of the isolates. All of the isolates lacking this mutation in clade B also have a second indel in the same locus 1,240 bases upstream.

Given the relatively large number of variants that distinguish clades A and B, it is impossible to conclusively identify the beneficial mutation that initially gave rise to clade B. Nevertheless, it may be possible to identify SNPs that helped maintain one or both populations. As discussed above, three independent nonsynonymous mutations arose during the course of this study in the pbpB locus, encoding PBP3. Two derived alleles are found exclusively in clade B isolates (one at very low frequency), and one derived allele is found exclusively in clade A isolates (Fig. 5; see Fig.  S5 and S6 in the supplemental material). The ancestral PBP3 allele is found in both clades.

FIG 5 

Distribution of pbpB alleles in clade A (A) and clade (B) over the 12 sampling time points. Clades A and B share an ancestral major allele. A minor allele is observed specifically in clade A, and two others segregate only in clade B.

We noted that aztreonam and piperacillin-tazobactam were first administered intravenously during the time that roughly coincides with the emergence of clade B (Fig. 1). Since aztreonam has high affinity for PBP3 (54), we examined the in vitro fitness of the clade A-derived pbpB allele, the most common clade B-derived pbpB allele, and the ancestral pbpB allele using standard MIC assays against four antipseudomonal β-lactam antibiotics (Fig. 6) (55). The ancestral allele is present in both clades but carried in different genetic backgrounds; nevertheless, the resistance it confers does not differ significantly among strains, regardless of their clade of origin (P = 0.15 to 0.44). The derived PBP3 alleles from both clades A and B increased strain resistance relative to the ancestral allele in the presence of both aztreonam and cefsulodin. In contrast, when tested on ceftazidime and piperacillin, while the clade A allele increased resistance, the clade B allele was no more resistant than the ancestral allele. These results indicate that the in vitro fitness contribution of each derived PBP3 allele (as measured by antibiotic resistance) is dependent on the specific antibiotic environment in which it is found.

FIG 6 

Relationship between pbpB alleles and MICs to selected antipseudomonal antibiotics. Pseudomonas aeruginosa isolates were grouped by pbpB genotype based on the presence of major and minor alleles, as well as clade membership in clade A (ancestral) and clade B (sweep) populations. Antibiotic MICs were assayed for aztreonam (A), cefsulodin (B), ceftazidime (C), and piperacillin (D). The presence of the minor alleles was associated with significantly increased resistance to aztreonam and cefsulodin in both clades but to ceftazidime and piperacillin only for the clade A minor allele. Statistical significance was determined by Mann-Whitney U test and is indicated with asterisks. (**, P ≤ 0.01; ***, P ≤ 0.001).

DISCUSSION

We have characterized the fine-scale evolutionary dynamics of P. aeruginosa during human CF lung infection and have observed a temporal shift in intraspecific population structure consistent with recurrent emergence and expansion of new clonal lineages. Since there is no change in the absolute size of the bacterial population over time, and Pseudomonas makes up the bulk of each of these populations (as shown by quantitative PCR [qPCR] analysis and microbiome sequencing, respectively), it is likely that the expansion of any new clonal lineage comes at the cost of preexisting clonal linages. Further, since there is strong evidence for positive selection in the population, we believe the best explanation for these data are recurrent selective sweeps, also known as periodic selection.

Although this patient has been infected with P. aeruginosa for at least 12 years, the entire population of P. aeruginosa coalesced to a common ancestor less than 3 years before the last specimen sampled for this study, indicating that this sweep went to fixation. It is possible that our sputum-based sampling did not capture the entirety of the P. aeruginosa population living within this patient and that there are additional distinct subpopulations of P. aeruginosa located in other anatomical compartments. However, if this were the case, then the migration rate between these subpopulations would have to be so low that we did not detect a single migrant over the year of sampling. Any rare migration from a divergent subpopulation would have been seen in the genealogical structure and time to most recent common ancestry of the overall population. On balance, we believe that such a low rate of migration between anatomical compartments is fairly unlikely and that we have therefore likely sampled the overall diversity of P. aeruginosa in this patient.

A second clonal lineage arose ~2 years after the first, giving rise to clade B. We cannot conclusively identify the “origin mutation” that gave rise to this lineage since there are a number of mutations that are fixed between the two clades. Consequently, no single variant is uniquely associated with the origin of clade B. This is perhaps not surprising since clade B originated in the nearly 6-month period between the collection of the first and second specimens. Any mutations that occurred in the emergent clade B common ancestor would be held in linkage disequilibrium with the causative mutation due to the very low rate of recombination.

Incomplete selective sweeps, such as those we observe in clade B, have been documented in other systems and generated substantial experimental and theoretical interest (56–69). A number of hypotheses can be proposed to explain why clade B has not swept to fixation. For example, we may have simply caught the sweep at an intermediate stage on its way to fixation, or the short-term persistence of the two clades could be due to stochastic processes. It is also possible that the ancestor of clade B colonized a distinct anatomical space and has since stayed isolated from its ancestral lineage. However, perhaps the most interesting explanation is clonal interference (70) driven by either competing beneficial mutations, negative epistasis between beneficial and deleterious mutations, or changing selection pressures driven by environmental heterogeneity. An example of the latter can be seen in our PBP3 data. While the clade A PBP3-derived allele confers a fitness advantage relative to the ancestral allele under all antibiotics tested, the clade B-derived allele is more fit than the ancestral allele only under half of the tested antibiotics. This is an example of latent potential for selection (71), in which the selective advantage of an allele is dependent on the environment. Even if the clade B allele is beneficial under certain conditions (such as during aztreonam treatment), it may be disadvantageous in other environments (such as during ceftazidime treatment), thereby, promoting the maintenance of the ancestral population and potentially explaining an underlying mechanism for clonal interference.

The analysis of natural selection within population data is complicated by short-term segregation of genetic polymorphisms, which can confound the more appropriate measure of fixed differences between populations (72). We have attempted to compensate for this issue by comparing dN/dS signals that emerge in the full data set versus the data set containing only the multispecimen SNPs. While the full data set provides a dN/dS value that is consistent with the expectation of a purely random mutational process (both synonymous and nonsynonymous mutations occurring at a rate consistent with the number of potential sites), the SNPs that are recovered from multiple specimens (multiple time points) over the course of the study show a signal consistent with strong positive selection in clade B, as would be expected in the case of a selective sweep. Again, it is not possible to use this analysis to identify the causative variant since the sweep population has had very little time for recombination to disentangle the selected substitution from the genomic background.

While dN/dS values can be difficult to interpret, a particularly notable indicator of selection is the number of loci that have accumulated multiple mutations over the short time since the most recent common ancestor arose. Our population of 233 strains has seven Mpl alleles, six AlgG and MucA/AlgU alleles, five AmpR/AmpC alleles, four PBP3 and MexA/MexB alleles, and three MexT and MucA alleles segregating (Table 3). The probability of even two mutations occurring at the same locus (giving rise to three alleles) is 1.64 × 10−4, or 0.027 when Bonferroni corrected for 165 independent mutations (107 SNPs and 58 indels). This pattern of multiple independent mutations is a powerful signal of parallel pathoadaptation (73) and mirrors the findings of Marvig and colleagues (28), who found extensive evidence for parallel evolution of P. aeruginosa in the CF lung. A number of multiply mutated loci overlapped between studies, including those that encode Mpl, MucA, AlgU, MexA/MexB, and LasR. A recent study by Williams et al. (29) also identified the mpl locus as carrying homoplasious (parallel) mutations in multiple strains, possibly due to homologous recombination. These patterns of parallel pathoadaptation are particularly striking given that nearly all of the affected loci are involved in either antimicrobial resistance or biofilm production.

The use of deep population sequencing in a single CF patient over a short period reveals community dynamics that are not evident by 16S rRNA gene microbiome sequencing or long-term longitudinal studies. Analyses of CF microbiota using culture-independent techniques (which are typically limited to genus-level taxonomic identification) have demonstrated that Pseudomonas commonly dominates microbial communities, is highly resilient, is associated with poor lung function and disease aggressiveness, and is a determinant of community dynamics (5, 6, 74–77). However, these analyses have been limited in their ability to predict temporal variation in patient status, such as the origins of exacerbations, and have not been able to explain the apparent stability of Pseudomonas through treatment. Our analysis shows that, although apparently stable at the genus level, Pseudomonas (P. aeruginosa specifically) can be a highly dynamic pathogen, rapidly adapting to environmental stress, which thus enables it to persist and resist the intended effects of antibiotic therapy.

The rapid evolution of P. aeruginosa revealed by deep sequencing in this study may help explain why significant changes in the antimicrobial treatment regimen did not result in any significant change in bacterial community composition. The apparent microbiologic stability observed belies a highly dynamic P. aeruginosa population that carries extensive genetic variation in a range of virulence- and resistance-associated loci, which are under strong positive selection. These population dynamics, which are undetectable using lower-resolution methods (e.g., microbiome- or culture-based methods), allow the pathogen to rapidly respond to therapeutic changes and could explain why some changes in treatment do not result in clinical improvement.

The PBP3 data reveal that changes in the selective environment (e.g., different antibiotics) can influence the fitness of strains in unpredictable ways. This observation may even be extended to the microbiome level, where changes in the relative abundance of Pseudomonas do not correlate with the anticipated response to antibiotics administered based on their antipseudomonal properties (such as the drop in Pseudomonas relative abundance following metronidazole administration at time point 9). Ultimately, we are left with a number of critical questions about the relationships between disease state, treatment approaches, and pathogen evolution. To what extent do irregular and variable antimicrobial treatments drive pathogen population dynamics? Are highly dynamic pathogen populations associated with specific clinical prognoses? Are there population structures or pathogen lineages that would favor better clinical outcomes? Is it possible to select for these lineages, and how could this be accomplished?

This study is limited to a single bacterial species in one infected individual over a defined time period, and thus its results may not generalize to other CF patients, disease stages, or bacterial species. Selective stress in human infection is highly heterogeneous and includes a broad array of differences in host response, microbial milieu, and treatment administration that would be expected to produce various temporal and spatial selective pressures, resulting in multiple patterns of intraspecies evolution. Future efforts to characterize the prevalence and dynamics of different patterns of selection across pathogens, individuals, and disease states are needed to understand their significance in health, disease, and treatment response.

MATERIALS AND METHODS

Specimen collection and microbiological analysis.

All protocols for the collection and use of patient sputum were approved by the Research Ethics Boards at the University Health Network (Toronto) and St. Michael’s Hospital (Toronto). Sputum samples were voluntarily produced by expectoration over a 350-day period by a single female study patient, CF67. Sputa were processed and P. aeruginosa cultured as previously described (40).

Following 72 h of incubation, all cultures were visually inspected and sampled for overall diversity in colony morphology. Twenty colonies resembling P. aeruginosa that were oxidase positive and exhibited growth on cetrimide and MacConkey agars (Becton, Dickinson, Hunt Valley, MD) were identified and cryopreserved for further analysis as previously described (38).

To study the bulk P. aeruginosa population, patient sputum was also spread plated directly onto P. aeruginosa-selective cetrimide agar and incubated at 37°C for 72 h. The resulting lawn of P. aeruginosa colonies was scraped from each replicate plate using a sterile inoculating loop and pooled in a single tube of sterile H2O for population sequencing.

DNA isolation.

Genomic DNA was extracted for whole-genome sequencing from 233 individual P. aeruginosa isolates and the pooled populations using the DNeasy blood and tissue kit (Qiagen, Ontario, Canada) as per the manufacturer’s protocol for Gram-negative bacteria.

Antimicrobial susceptibility testing of pbpB variants.

MICs for the β-lactam antibiotics aztreonam (Alfa-Aesar, Ward Hill, MA), ceftazidime (Sigma-Aldrich, Ontario, Canada), cefsulodin (A.G. Scientific, Inc., San Diego, CA), and piperacillin (Sigma-Aldrich, Ontario, Canada) were determined for 60 pbpB variants (20 per genotype) by agar dilution, as described by the Clinical and Laboratory Standards Institute guidelines (78). Significant differences in fitness defined by variations in MIC associated with clade A or clade B isolates carrying the pbpB ancestral major allele, ancestral minor allele, or sweep minor allele were assessed in GraphPad Prism 5.0 (GraphPad Software, Inc., San Diego, CA) using a Mann-Whitney U test with a 95% confidence interval.

Microbiome analysis and 16S quantitation.

Microbiome analysis was performed by analysis of the V5 to V7 hypervariable regions of the 16S rRNA gene using the SI-Seq protocol (38). Sequencing was performed at the University of Toronto, Centre for the Analysis of Genome Evolution and Function (CAGEF). An in-house Java script was used to prepare the sequence for downstream analysis using MACQIIME (79). Chimeric reads were removed using reference-based chimera detection with USEARCH (80). Operational taxonomic units (OTU) were picked against an SI-Seq structured reference set using a similarity of 0.87 (corresponding to the genus level by the SI-Seq method), with premature termination turned off. Taxonomy was assigned using an RDP structured data set using UCLUST with --max_accepts set to 0 and similarity set to 0.97 (empirically determined to produce the highest consistency identification at the genus level using P. aeruginosa controls) (80, 81). Unassigned OTU were identified with BLASTN (82). OTU representing less than 0.005% of overall sequence were removed. Quantitative PCR of the 16S rRNA gene was performed as previously described (83).

Sequencing and quality control.

Whole-genome sequencing of sampled P. aeruginosa isolates was performed on the GAIIx, HiSeq, and MiSeq Illumina platforms. The distribution of the number of bases sequenced per isolate ranged from 167 to 1,426 million bases, and the median was 402 million bases (see Fig. S1A in the supplemental material). This sequencing depth yielded a median read depth per position of 54× (range 25× to 183× [see Fig. S1B]). Sequencing reads were quality trimmed employing Quake and its native k-mer counter (84), and adapters were removed in the CLC Genomics Workbench v. 6.5.1 (Aarhus, Denmark).

De novo and reference mapping assembly.

Isolate CF67-71, which was sequenced with 250-bp-long paired-end reads, was assembled de novo using the CLC Genomics Workbench (Aarhus, Denmark), using its default assembly settings. Contigs with a scaffolding depth of lower than 10× and/or with a size smaller than 1 kb were removed from this analysis. Contigs that passed this filter were annotated at the RAST server using the native gene caller and Classic RAST as the annotation scheme (85). CF67-71 contigs were used as the reference for mapping assembly of each remaining isolate. We performed six different reference-mapping assemblies using BWA v. 0.7.4 (86), LAST v. 284 (87), and novoalign v. 2/8/03 (Novocraft Technologies), using sequencing reads with and without quality trimming (Fig. 2).

SNP and indel calling.

The preliminary list of variants was extracted from the reference-mapping step using SAMtools and BCFtools v. 0.1.19 (88). We then employed a two-step SNP calling process similar to that of Lieberman et al. (89). In the first step, we generated a set of 108 polymorphic positions by identifying the overlapping set of variants called by all six reference-mapping approaches using a custom Java pipeline (available at https://github.com/DSGlab/SNPCallingPipeline). This step used the following criteria: (i) variant Phred quality score of ≥30 and (ii) variants must be found at least 150 bp away from the edge of the reference contig and/or an indel. In the second step, we revised each polymorphic position from the previous step in all of the isolates and included calls that had a Phred score of at least 25. Finally, we tested the sequencing support for either the SNP or the reference base. This multihypothesis correction required that at least 80% of the reads supported the SNP or the reference. If the data did not support either base, then the position was called as an ambiguous base (“N”).

We used Dindel v. 1.01 (90) to realign the sequencing reads (mapped and unmapped) and test whether candidate indels produced by BWA and SAMtools were supported by those reads. Putative indel positions were selected as having an indel call in at least one nonreference isolate with the following criteria: (i) a variant Phred quality score of ≥35, (ii) at least 2 forward and 2 reverse reads, and (iii) sequencing coverage of ≥10. Each of these putative indel positions was then reviewed in each of the 233 isolates. An indel was called if the variant had a Phred quality score of ≥25 and the allele frequency was ≥80%. Indel calls where the allele frequency was ≤20% were labeled as ambiguous. All other indel calls were discarded.

Population and single-genome-sequencing evaluation.

To confirm that our sampling design adequately sampled the underlying P. aeruginosa population, we compared the major allele frequency for each segregating site from individual isolates with that obtained from bulk population sequencing of P. aeruginosa isolates collected from 1 of the 12 specimens (see Fig. S3 in the supplemental material).

Sequencing reads from each isolate from specimen 9 were rarified to 1/20th of the total number of reads in the population sequencing experiment such that the overall number of reads would be the same. Reads from both experiments were mapped to the contigs of isolate CF67-71 as described above. Major and minor allele frequencies were calculated using the quality thresholds described by Lieberman et al. (89).

Phylogenetic, population genetic, and coalescent analyses.

An alignment of the 107 SNPs was employed to reconstruct the phylogeny of the 233 isolates by both maximum likelihood and neighbor-joining methods using maximum composite likelihood distances and gamma correction (α = 1) in MEGA6 (91). The two trees gave nearly identical results. The phylogeny was also investigated with a network-based phylogenetic analysis. SplitsTree4 (92, 93) with the Jukes Cantor distance matrix was used to create a neighbor-net phylogenetic tree.

The time to the most recent common ancestor (tMRCA) among our 233 isolates and between our isolates and the P. aeruginosa strain LESB58 was calculated by Bayesian Markov chain Monte Carlo (MCMC) analysis using BEAUTi and BEAST v. 1.75 (43, 44). The parameters deduced from prior information were days between isolation dates and nucleotide substitution model and frequencies. The nucleotide substitution model that best fit our data according to FindModel (http://www.hiv.lanl.gov/content/sequence/findmodel/findmodel.html), which is a web application of ModelTest (94), was the general time-reversible (GTR) model (Log likelihood [LnL] = −927.007482, and Akaike information criterion [AIC] = 1,870.014964). The nucleotide substitution frequencies AC = 0.0661, AG = 0.3306, AT = 0.0275, CG = 0.1858, GT = 0.0657, and CT = 0.3241 were estimated employing the Maximum Likelihood Estimate of Substitution Matrix tool in MEGA6 (91). Preliminary MCMC analyses consisting of 10 million generations were employed to infer the appropriate molecular clock and demographic models (95). We inferred the log-normal relaxed uncorrelated molecular clock for the constant size coalescent analysis. Given the estimated priors, the final analysis was run in duplicate through 1 billion MCMC generations sampled every 1,000 MCMC generations, and the burn-in period was set at 100 million MCMC generations. The inferred molecular clock was consistent with the number of mutations observed in our isolates through time (see Fig. S7 in the supplemental material).

Population genetic tests were performed with DnaSP v. 5.10.01 (96). P values were determined via coalescent simulations as implemented in DnaSP.

Interstrain whole-genome alignment and phylogeny.

The full genomic sequences of the P. aeruginosa strains LESB58, PAO1, DK2, PA7, NCGM2, UCBPP, B136, and M18 were obtained from the Pseudomonas Genome Database (97). These sequences and the contig sequences of the de novo-assembled isolate CF67-71 were aligned using Mugsy v. 1.2.3 (98). The output alignment file in MAF format was converted to the more broadly used FASTA format in the Galaxy platform (99). Finally, the neighbor-joining phylogeny (see Fig. S4 in the supplemental material) was inferred using the CLC Genomics Workbench v. 6.5.1 (Aarhus, Denmark) with a bootstrap value of 500.

Whole-genome sequencing of 233 isolates of P. aeruginosa with Illumina. (A) Distribution of number of bases sequenced per isolate. The distribution ranges from 167 to 1,426 million sequenced bases, and the median is 402 million sequenced bases. (B) Distribution of median read depth per position. Download

Figure S1, TIF file, 2 MB

Flowchart of the SNP calling pipeline. All sequencing reads were quality trimmed and assembled de novo. The isolate with the highest N50 was used as the reference to which the sequencing reads of all isolates were mapped using three algorithmically distinct methods. SNPs called by all of these methods were selected to avoid biases from sequencing quality or coverage. Download

Figure S2, TIF file, 2.6 MB

Individual isolate sequencing reflects whole-population sequencing adequately. The frequency of the major allele of 21 SNPs among the 20 isolates collected in the 9th time point is plotted against the frequency of the alleles as determined by population sequencing of the specimen. No alleles were recovered in population sequencing that were not found in the individual isolate sequencing. Download

Figure S3, TIF file, 1.5 MB

CF67 isolates in the context of other P. aeruginosa complete genomes. Shown is the unrooted neighbor-joining (NJ) phylogeny of CF67 isolate 71 and 8 P. aeruginosa complete genomes based on genome-wide alignment. All nodes have 100% bootstrap support. Download

Figure S4, TIF file, 1.2 MB

Maximum likelihood tree annotated with PBP3 alleles. The tree is constructed as described in the legend to Fig. 2. The identity of the PBP3 allele carried by each strain is indicated by the color of the box to the left of the isolate name. Download

Figure S5, TIF file, 1.4 MB

Location of nonsynonymous mutations in the pbpB gene in the 3D crystal structure of PBP3. (A) Multiple sequence alignment of the PBP3 protein sequences of the ancestral allele, minor allele of clade A, and the most common minor allele of clade B. The mutated amino acid is shaded in red and green for clades A and B, respectively. (B) The crystal structure of PBP3 complexed with aztreonam (PDB ID 3PBS [42, 43]) was customized on the PyMOL Molecular Graphics System v. 1.5.0.4 (Schrödinger, LLC). The aztreonam structure is in white, and the mutated positions are colored in red for the clade A minor allele and green for the clade B minor allele. Download

Figure S6, TIF file, 1.7 MB

Temporal signal compared to the inferred substitution rate of the 233 isolates. (A) Regression analysis of the root-to-tip distance as a function of time of isolation using the Path-O-Gen program (http://tree.bio.ed.ac.uk/software/pathogen/) as described in reference 101. Each circle represents the average root-to-tip distance of the isolates from the respective sampling time point. The resulting trend shows that the inferred molecular clock was consistent with the changes seen in our isolates through time (R2 = 0.79, P = 0.0001). (B) Rate of substitution as inferred by the coalescent analysis (2.92 ± 0.70 SNPs per year). The line is set to intersect the x axis at the tMRCA of all isolates of this study. Download

Figure S7, TIF file, 1.2 MB

Patient details for the study period.

Table S1, DOCX file, 0.01 MB

SNPs implicated in antibiotic resistance and host adaptation.

Table S2, DOCX file, 0.02 MB