ABSTRACT
Pseudomonas aeruginosa is an opportunistic human pathogen that frequently causes health care-associated infections (HAIs). Due to its metabolic diversity and ability to form biofilms, this Gram-negative nonfermenting bacterium can persist in the health care environment, which can lead to prolonged HAI outbreaks. We describe the creation of a core genome multilocus sequence typing (cgMLST) scheme to provide a stable platform for the rapid comparison of P. aeruginosa isolates using whole-genome sequencing (WGS) data. We used a diverse set of 58 complete P. aeruginosa genomes to curate a set of 4,440 core genes found in each isolate, representing ∼64% of the average genome size. We then expanded the alleles for each gene using 1,991 contig-level genome sequences. The scheme was used to analyze genomes from four historical HAI outbreaks to compare the phylogenies generated using cgMLST to those of other means (traditional MLST, pulsed-field gel electrophoresis [PFGE], and single-nucleotide variant [SNV] analysis). The cgMLST scheme provides sufficient resolution for analyzing individual outbreaks, as well as the stability for comparisons across a variety of isolates encountered in surveillance studies, making it a valuable tool for the rapid analysis of P. aeruginosa genomes.
INTRODUCTION
Pseudomonas aeruginosa is a ubiquitous Gram-negative bacterium responsible for 32,600 multidrug-resistant, health care-associated infections (HAIs) per year in the United States (1). P. aeruginosa can cause a variety of HAIs, including pneumonia and bloodstream and wound infections; it is also the main cause of mortality in cystic fibrosis patients (2).
P. aeruginosa has a large, complex genome with intrinsic, adaptive, and acquired resistance mechanisms (3) and multiple virulence factors (4), making infections both difficult to treat and highly morbid (5). P. aeruginosa can also form biofilms, allowing it to persist in the health care environment, including on surfaces (6), medical devices (7), and in water sources (8). Because of its ubiquity, outbreaks involving P. aeruginosa can be both clonal and involve multiple sequence types (STs) (9), as defined by the traditional seven-gene multilocus sequence typing (MLST) scheme (10, 11).
The advent of whole-genome sequencing (WGS) has changed HAI outbreak investigations by drastically increasing the resolution with which isolates can be genetically characterized (12, 13). Traditional MLST schemes can be expanded from seven genes to core genome MLST (cgMLST) schemes with WGS data by including the thousands of genes common to a particular species (14) or on an ad hoc basis (15). Herein, we describe the creation of a cgMLST scheme for P. aeruginosa and apply it to explore the diversity of the species in publicly available genomes. We also compare it to other phylogenetic methods using isolates from a convenient set of four HAI outbreaks.
MATERIALS AND METHODS
Development of the P. aeruginosa cgMLST scheme.Genes from all 58 complete P. aeruginosa genomes available from RefSeq (16) in January 2017 were aligned and compared using BLAT v36x2 (17). All genes present as a single copy in each genome with ≥95% nucleotide similarity to the allele from the reference PAO1 genome were used to construct a core genome MLST (cgMLST) scheme. The cgMLST scheme was implemented using custom scripts (https://github.com/DHQP) on an additional 1,991 genomes downloaded from RefSeq (16) for evaluation, to refine the loci included in the scheme, and to expand the alleles for each gene. A Bonferroni correction threshold of P < 0.05 (P = 0.05/2,049 = 1.1 × 10−5) was used to determine whether any genes should be removed from the cgMLST scheme if they were missing, present as multiple copies, or varied by ≥10% in size at a statistically significant rate in the total set of 2,049 genomes (58 complete plus 1,991 assemblies). This was implemented by determining a Z-score associated with the absence of each gene, and genes with Z-scores higher than that of P were removed. For each gene in the scheme the number missing, present as multiple copies, and Z-score for all RefSeq genomes analyzed are included in Spreadsheet S2 in the supplemental material.
Identification of prophage and recombination genes.The web server PHASTER (18) (phaster.ca) was used to identify genes included in the cgMLST scheme that originated form prophage sequences (Spreadsheet S2). ClonalFrameML (19) was run using default parameters on a dendrogram made with the cgMLST scheme for the 58 complete P. aeruginosa genomes to find genes that included regions of recombination (Spreadsheet S2).
Tree generation from cgMLST scheme.DendroPy 4.2 (20) was used to create dendrograms using the unweighted pair group method with arithmetic mean (UPGMA) and neighbor joining (NJ) methods based on the number of allelic differences in the cgMLST scheme between genomes. The UPGMA method generated dendrograms more similar to the trees generated by the hqSNV pipeline (Table S2), so was used as the primary method in this work. Bootstrap support values were determined using the SumTrees package in DendroPy.
Pulsed-field gel electrophoresis.Pulsed-field gel electrophoresis (PFGE) patterns were determined using standard procedures for Gram-negative bacteria (https://www.cdc.gov/hai/pdfs/labsettings/Modified-PulsedNet-procedure-GNB.pdf) and analyzed through the BioNumerics software package from Applied Maths NV (Sint-Martens-Latem, Belgium). DendroPy 4.2 (20) was used to create dendrograms using UPGMA clustering based on the number of band differences in the PFGE patterns between isolates.
Whole-genome sequencing.DNA was extracted using the Promega Maxwell 16 Cell Low Elution Volume (LEV) DNA purification kit and the automated Maxwell 16 MDx instrument (Madison, WI). High quality input genomic DNA (gDNA) was fragmented with the Covaris ME220 Focused-ultrasonicator (Woburn, MA). Sample libraries were prepared using the NuGEN Ovation Ultralow System V2 assay kit (San Carlos, CA). Sequencing was completed using the Illumina MiSeq platform (San Diego, CA) to produce 250-bp paired-end reads.
Genome assembly, annotation, and sequence typing.Genomes were assembled de novo using high quality reads with SPAdes v3.9 (21). Genes were identified using Prokka v1.12 (22). Sequence types (STs) were determined using the multilocus sequence typing (MLST) definitions from PubMLST (https://pubmlst.org/paeruginosa/).
Generation of SNV trees.The hqSNV pipeline SNVPhyl 1.3.0 (23) was used to generate phylogenetic trees using sequenced reads from the isolates associated with the four outbreaks, with the SNV abundance set to 0.75, minimum coverage set to 10, and the filter-density threshold and window set to 2 and 11, respectively. The median assembly based on Mash distance (24) was used as the reference. The k-mer-based trees were generated with the SNV pipeline kSNP v3.0 (25), using the assembled genomes from the outbreaks and a k-mer length of 13.
Tree comparisons.Statistical comparisons of PFGE, cgMLST, and SNV tree topologies generated from the outbreak isolates were made using the unweighted Kendall-Colijn similarity metric implemented through the R package treespace (26). P values were calculated for the similarity of trees using a Z-test on a population standard deviation generated from 10,000 comparison scores from random trees (27). A threshold of P < 0.05 was used to determine similarity. Tree images were created and annotated using iTOL v3 (28).
Data availability.Raw sequences have been deposited in NCBI BioProject ID PRJNA288601.
RESULTS
Development of the P. aeruginosa cgMLST scheme.All 58 complete P. aeruginosa genomes available in NCBI’s RefSeq (16) in January of 2017 were used to construct the initial cgMLST gene set, as depicted in Fig. 1. There were 36 unique STs represented in the set, which also included eight genomes from high-risk clones (29), composed of three each of ST111 and ST253 and two from ST235 (Spreadsheet S1). The average genome size was 6.7 MB with 6,096 open reading frames. Three of the sequences included plasmids, which were excluded from the construction of the cgMLST scheme. There were 353,592 total gene alleles in the entire set, which when aligned represented 17,258 gene loci; 4,514 of the genes were present as a single copy in each genome and were included in the initial cgMLST gene list, which was later reduced to 4,440 genes after assessing an additional 1,991 contig-level genomes. These genes totaled 4.3 MB, representing ∼64% of the average P. aeruginosa genome.
Creation of P. aeruginosa cgMLST gene set.
The core genome included 2,780 genes with known or probable functions, including the seven from the traditional MLST scheme. There were five virulence factors (exoT, lasB, lecA, toxA, and toxR) (4), as well as 23 genes involved in antimicrobial resistance, including one fosfomycin resistance enzyme, two beta-lactamases, and 20 components of multigene, multidrug efflux pump systems (30) (Spreadsheet S2).
The core gene list was also compared to representative genomes from the other 12 members of the P. aeruginosa species group (NCBI:txid136841) (31). They contained an average of 148 (3.4%) of the cgMLST genes, where Pseudomonas citronellolis had the most with 463 (10.5%), while Pseudomonas caeni had no overlapping genes (Table S1).
Exploration of P. aeruginosa diversity.The 1,991 contig-level P. aeruginosa genome assemblies available in RefSeq (16) as of January of 2017 were assessed using the cgMLST scheme to explore the diversity of the species. These included 336 unique STs, 154 with more than one representative genome, as well as 190 genomes without defined STs.
Because assembled genomes may have artifacts of short-read sequencing, such as missing, incomplete, or duplicated genes (13), it was expected that some of the cgMLST loci would be absent or overrepresented from isolates in the expanded set. One hundred and fourteen genes were removed from the scheme because they were either missing from or present as multiple copies in too many genomes (see details in the Materials and Methods section), leaving a final set of 4,440 loci (Spreadsheet S2). On average there were 21 genes missing per sample, giving a coverage of 99.5% of the 4,440-gene cgMLST scheme, and in the average assembly there were 1.4 alleles with multiple copies (Spreadsheet S3).
We examined the core gene set for evidence of homologous recombination by applying ClonalFrameML (32) to a dendrogram made from the 58 complete genomes. There were 510 unique recombination sequences identified on 369 of the genes, representing 0.55% of the 92,837 unique cgMLST alleles in those genomes. Two-thirds (342/510) of the unique recombination sequences were found on multiple genomes, supporting the observation that homologous recombination helps shape the core genome of bacteria (33). Since recombination has been to shown to have greater effect on branch length rather than tree topology (34), which is of more concern when assessing relatedness in HAI outbreaks, and the alleles affected still met the statistical criteria for inclusion (Spreadsheet S2), they were not removed from the scheme.
Additionally, 20 (0.45%) of the 4,440 core genes were similar to prophage sequences, as determined by PHASTER (18). While prophage sequences can have higher rates of mutation, they can also accelerate evolution when stably incorporated into the genome (35), and thus contribute to phylogeny. Like the recombination alleles, there was no empirical evidence of increased variation in the prophage alleles in the cgMLST scheme (Spreadsheet S2), so they too were kept in the scheme.
A pairwise comparison between the alleles found in each isolate was used to make a difference matrix for all isolates, from which the UPGMA was used to create a dendrogram (Fig. 2). The isolates generally clustered by ST, as different STs varied by an average of 3,889 alleles (median: 3,739; range: 8 to 4,300), while genomes from the same ST differed by an average of 331 alleles (median: 209; range: 0 to 2,903). Nineteen pairs of genomes differed by less than 10 alleles, including four pairs with zero allele differences. However, there were also instances of large diversity within STs, as 27 included isolate pairs that differed by more than 1,000 alleles. And the STs were not completely distinct in all cases, as 47 genomes had more common alleles with genomes from different STs rather than their own. Of course, some of the extreme observed differences within STs could be due to variations in sequence quality.
Dendrogram from cgMLST scheme of 2,049 P. aeruginosa genomes. Stars indicate entries used to construct the core gene list. The most common MLST types are highlighted on the dendrogram.
Analysis of isolates from HAI outbreaks.Isolates from four HAI outbreaks were used to compare the phylogeny created using the cgMLST scheme to those made from pulsed-field gel electrophoresis (PFGE) patterns and single nucleotide variant (SNV) pipelines. Outbreak 1 occurred in a single hospital unit (36) and included 12 isolates from three STs; outbreak 2 occurred in a skilled nursing facility with ventilated residents (9) and included 28 isolates from four different STs; outbreak 3 occurred in an intensive care unit and included nine isolates from three different STs; and outbreak 4 included 26 isolates from a single ST that were found across multiple facilities within the same geographic area. Only clinical isolates were included in this analysis.
PFGE, a means for determining relatedness between bacterial genomes using restriction digestion (37), generates dozens of fragments and offers more discriminatory power than traditional MLST for P. aeruginosa (38). Isolates from outbreaks 1, 2, and 3, for which both PFGE and WGS data were available, were selected for comparison. Only two isolates from outbreak 4 underwent both PFGE and WGS, so it was excluded from this analysis. Differences in PFGE patterns were used to construct UPGMA dendrograms. These were compared to trees created by the cgMLST scheme using the Kendall-Colijn metric (26), a quantitative measure of the difference in topology between two phylogenetic trees (27).
The topologies of the PFGE and cgMLST trees were statistically similar for outbreaks 1, 2, and 3 (Table S3). The tree topologies also showed similar clustering by ST for outbreaks 1 and 3, and by facility in outbreak 2 (as all the isolates were from the same ST) (Fig. 3), validating the results obtained from the statistical comparison. The isolates with identical PFGE patterns varied by a median value of 26 (0.59%) and 29 (0.66%) alleles for outbreaks 1 and 3, respectively. There were no isolates with identical PFGE patterns from outbreak 2.
PFGE and cgMLST trees of isolates from three P. aeruginosa HAI outbreaks. Isolates from outbreaks 1 and 3 are colored by MLST sequence types (ST), while those from outbreak 2 (which were all from the same ST) are colored by the facility they were recovered from.
Given that P. aeruginosa genomes from the same ST vary by only a few hundred alleles on average, the cgMLST scheme was tested to see if it could differentially cluster outbreak isolates from genomes of the same ST from RefSeq. The median allele difference between the unrelated RefSeq and the outbreak isolates was more than four times greater than within the outbreak isolates themselves for all outbreaks (Table 1), and trees generated using the cgMLST scheme show distinct clusters separating the outbreak from RefSeq genomes (Fig. S1). There were no instances where an outbreak isolate was more closely related to a RefSeq genome of the same ST than to another outbreak isolate.
Differentiation of outbreak isolates versus unrelated RefSeq isolates of same sequence type
SNV-based WGS comparison techniques are commonly used in HAI outbreak investigations to create phylogenetic trees to determine isolate relatedness and possible sources of infection (39). While cgMLST can measure thousands of allele-level differences, SNV-based approaches consider millions of possible nucleotide-level variations in both coding and noncoding sequences, providing much greater resolution, though analysis can be complicated by factors such as recombination (19) and convergent evolution (40). Trees made from the cgMLST scheme were compared to those from two SNV tools: SNVPhyl, an hqSNV pipeline (23), and kSNP, a k-mer based SNV program (25). SNV applications rely on large aligned core genomes for comparison (41), so only closely related isolates from the four outbreaks were included in the analysis. For outbreaks 1, 2 and 4, isolates from the predominant ST were included. For outbreak 3, the cgMLST tree showed that the isolate from ST2775 was similar to the ST309 isolates (Fig. 3), so it was also included in the SNV analysis.
The dendrograms made from the cgMLST distances had quantitatively similar topologies to the trees from the k-mer based SNV pipeline for outbreaks 1, 3 and 4, and for outbreaks 3 and 4 for the hqSNV pipeline (Table 2). The trees from outbreak 4 showed similar clustering, as illustrated in the same two subgroups emerging in both the cgMLST and hqSNV trees (Fig. 4). The cgMLST allele differences and SNV counts correlated with each other as well, though there was much more variability for outbreaks 2 and 4, which both had more than twice as many isolates as the other two outbreaks (Fig. S2). Additionally, there were drastic differences in the slope of the correlations between SNVs and cgMLST allele differences for the k-mer and hqSNV pipelines, because areas of homologous recombination were filtered out by the hqSNV pipeline but included in the cgMLST scheme (Table S4) and the k-mer pipeline. While these regions affected only a fraction of the unique alleles in the outbreak isolates (0.76 to 3.8%), they are overrepresented in the k-mer pipeline SNV counts because of the higher density of SNVs in these alleles resulting from recombination.
Comparison of cgMLST to SNV Trees
Trees generated by cgMLST scheme and hqSNV pipeline (SNVPhyl) for outbreak 4 isolates. Both techniques revealed two clusters of the same genomes (outlined in red and blue) distinct from the larger set. Support values for 1,000 bootstraps are labeled on the cgMLST tree.
DISCUSSION
Because its pathogenicity is not limited to specific strains (42) and multiple STs are often found in the health care environment (9), investigations of outbreaks of P. aeruginosa HAIs require flexible classification tools that can rapidly assess relatedness among both diverse and genetically similar isolates. The P. aeruginosa cgMLST scheme presented herein provides a means for comparison of WGS data using a simple, static scheme that can easily accommodate diverse samples and, as well, rapidly add new isolates into the analyses.
As expected, the cgMLST scheme generated clusters consistent with traditional, seven-gene MLST typing, but the increased resolution from the thousands of genes included in the scheme revealed heterogeneity within STs when applied to thousands of publicly available P. aeruginosa genomes. The cgMLST scheme revealed 47 genomes in the public set with more similarity to members of different STs than their own. We also observed this same phenomenon in the isolates from outbreak 3, as a group of four ST309 isolates had more common alleles with an isolate from ST2775 than with the other ST309 isolates from the same facility. While ST2775 and ST309 are closely related (differing only by one of seven MLST alleles), an investigation of this outbreak using MLST as a filter of relatedness would exclude the ST2775 isolate.
PFGE, another tool traditionally used in outbreak investigations, incorporates dozens of fragments to provide more resolution than MLST, but still much less than cgMLST schemes, and simplifies phylogenetic relationships versus sequencing-based approaches. This was evident when comparing PFGE-derived dendrograms to those made using the cgMLST scheme for three P. aeruginosa HAI outbreaks. While both techniques generated similar clustering, the cgMLST trees revealed subtle relationships between the isolates, as expected from the increase in resolution.
While both approaches utilize WGS data, SNV-based techniques offer a more granular means of genomic comparison than cgMLST. But because they abstract thousands of gene sequences to simple allele names, cgMLST schemes are less computationally demanding than the higher resolution SNV-based techniques, and cgMLST analyses do not have to be rerun on the entire data set when new isolates are added because the core genome is static, making them ideal for longitudinal analysis of large collections (43–45). When applied to the outbreak isolates, the scheme produced topologically similar phylogenetic trees compared to an hqSNV pipeline for two of the four outbreaks. This variability is likely due to multiple factors in addition to the fundamental difference of comparing SNVs to cgMLST alleles, such as the larger aligned core genome size generated by SNV pipelines (>93% of the genome for the hqSNV pipeline, which includes noncoding and accessory genomic elements) compared to the cgMLST scheme (∼64%), which only considers the coding regions of species-wide core genes. The hqSNV pipeline also explicitly filtered regions of recombination, some of which were included in the cgMLST scheme. However, recombination events affected only a subset of core genes and, in general, would not have outsized influence on tree topology since the scheme considers any change to a single gene to be equivalent, whether it arises from one, ten, or hundreds of SNVs. The method of tree construction also differs between approaches (maximum likelihood for the hqSNV pipeline and UPGMA for the cgMLST scheme). These differences are largely a result of the cgMLST scheme’s simplicity, but this same feature allows for increased speed of execution, making it a useful first pass filter to identify outbreak clusters that may require higher resolution SNV analysis.
The cgMLST scheme for P. aeruginosa provides a standardized means of rapid typing and comparison for this diverse, ubiquitous pathogen. It can be applied to surveillance efforts to assess the population structure of the pathogen on larger scales and help uncover previously undetected clusters of related infections before outbreaks occur in a facility or region. In addition, because P. aeruginosa health care-associated outbreaks often involve unrelated and nonclonal samples, it can rapidly identify potential transmission events for further analysis. Given these applications, the P. aeruginosa cgMLST scheme is a useful addition to the public health toolkit for WGS analysis. The scheme has been adapted by Applied Maths NV (Sint-Martens-Latem, Belgium) for use on their BioNumerics software package.
ACKNOWLEDGMENTS
The findings and conclusions in this report are those of the authors and do not necessarily represent the views of the Centers for Disease Control and Prevention. The use of trade names is for identification only and does not imply endorsement by the Centers for Disease Control and Prevention.
FOOTNOTES
- Received 4 February 2020.
- Returned for modification 2 March 2020.
- Accepted 28 May 2020.
- Accepted manuscript posted online 3 June 2020.
Supplemental material is available online only.
- This is a work of the U.S. Government and is not subject to copyright protection in the United States. Foreign copyrights may apply.
