Comparative Genomics of Nonoutbreak Pseudomonas aeruginosa Strains Underlines Genome Plasticity and Geographic Relatedness of the Global Clone ST235.

Pseudomonas aeruginosa is an important opportunistic pathogen in hospitals, responsible for various infections that are difficult to treat due to intrinsic and acquired antibiotic resistance. Here, 20 epidemiologically unrelated strains isolated from patients in a general hospital over a time period of two decades were analyzed using whole genome sequencing. The genomes were compared in order to assess the presence of a predominant clone or sequence type (ST). No clonal structure was identified, but core genome-based single nucleotide polymorphism (SNP) analysis distinguished two major, previously identified phylogenetic groups. Interestingly, most of the older strains isolated between 1994 and 1998 harbored exoU, encoding a cytotoxic phospholipase. In contrast, most strains isolated between 2011 and 2016 were exoU-negative and phylogenetically very distinct from the older strains, suggesting a population shift of nosocomial P. aeruginosa over time. Three out of 20 strains were ST235 strains, a global high-risk clonal lineage; these carried several additional resistance determinants including aac(6')Ib-cr encoding an aminoglycoside N-acetyltransferase that confers resistance to fluoroquinolones. Core genome comparison with ST235 strains from other parts of the world showed that the three strains clustered together with other Brazilian/Argentinean isolates. Despite this regional relatedness, the individuality of each of the three ST235 strains was revealed by core genome-based SNPs and the presence of genomic islands in the accessory genome. Similarly, strain-specific characteristics were detected for the remaining strains, indicative of individual evolutionary histories and elevated genome plasticity.

Introduction

Pseudomonas aeruginosa is a ubiquitous Gram-negative bacterium and one of the most predominant organisms associated with hospital-acquired infections (Vincent et al. 2009). Treatment of infected patients is further complicated due to the intrinsic resistance of this bacterium to therapeutics and its elaborated capacity to develop resistance to antibiotics through the selection of mutations in chromosomal genes and through horizontal gene transfer (Kulasekara et al. 2006; Kitao et al. 2012; Hong et al. 2015). P. aeruginosa response to environmental stress as well as its extensive metabolic capacity have been associated with the ability of the bacterium to persist in clinical settings (Frimmersdorf et al. 2010; Klockgether et al. 2011).

Even though P. aeruginosa has a nonclonal population structure, multilocus sequence typing (MLST) has revealed that a few sequence types (STs) are widely distributed and have been called high-risk clones. ST235 is the most important clonal ST lineage associated with outbreaks and poor clinical outcomes (Hong et al. 2016; Miyoshi-Akiyama et al. 2017; Treepong et al. 2018). Some important features of ST235 are resistance to fluoroquinolones, aminoglycosides and beta-lactams and the presence of virulence factors such as a type III secretion system (Finck-Barbançon et al., 1997; Sato et al. 2003; Kulasekara et al. 2006). Additionally, through whole-genome sequencing (WGS) studies in recent years, detailed insights into the population structure of P. aeruginosa have been reported. It is now evident that the population consists of three phylogenetically distinct clades or groups (Stewart et al. 2014; Freschi et al. 2015). Most P. aeruginosa strains sequenced to date belong to group 1, including the widely used reference strain PAO1 (Stover et al. 2000). Group 2 strains include the ST235 clonal linage, and are exoU-positive, in contrast to group 1 strains (Lee et al. 2006). Only few strains belong to group 3, including the reference strain PA7 (Roy et al. 2010).

The molecular epidemiology, resistance, and virulence properties of P. aeruginosa in nonoutbreak settings have received limited attention (Dettman et al. 2013; Hu et al. 2017). Here, we sequenced and analyzed epidemiologically unrelated P. aeruginosa strains, isolated from patients from a general hospital over a time period of two decades. We wanted to investigate to which degree these strains were related, if they belonged to a specific population group, clade or ST, and if detailed genome features of strains isolated at different times could indicate an environment pressure. Moreover, a comparison with global genomes was carried out in order to evaluate worldwide spread and molecular diversity of related clones.

Materials and Methods

The Clinical Laboratory from Hospital Israelita Albert Einstein, a Brazilian general hospital located in Sao Paulo, banks all isolates from normally sterile sites and all isolates obtained from documented nosocomial infections. Isolates are routinely identified to the species level using the Vitek 2 card system (Biomerieux, France). From this depository, 20 epidemiologically unrelated P. aeruginosa isolates collected over a period of 22 years were randomly selected for this study (table 1). Bacterial isolates were cultivated on Luria–Bertani broth for 18 h and genomic DNA was extracted using the DNeasy kit (Qiagen). Nextera XT shotgun libraries were prepared for the MiSeq (Illumina, San Diego, CA, USA) with a 300 bp paired-end sequencing run. Libraries were prepared according to the manufacturer’s protocol at the Göttingen Genomics Laboratory, Germany. Raw reads were quality controlled with FastQC v0.11.2 (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/; last accessed June 19, 2018), and subsequently trimmed using Trimmomatic 0.36 (http://www.usadellab.org/cms/? page=trimmomatic; last accessed June 19, 2018) to remove sequences with quality scores lower than 20 (Illumina 1.9 encoding) and remaining adaptor sequences. De novo assembly was performed using SPAdes v3.11 (Bankevich et al. 2012). Due to the Illumina sequencing technology used in this study, assembly problems related with repetitive sequences, including ribosomal RNA operons, transposases (often annotated as “mobile element protein”), and phage regions, are unavoidable. Low sequence read coverage at random positions is also often reported and may lead to assembly issues. We did not further investigate any particular gap or contig ends since our study focused on phylogenomic analyses. The ST of all strains analyzed in this work was determined from their whole genome sequence data using the multilocus sequencing typing (MLST) scheme (https://pubmlst.org/paeruginosa/; last accessed June 19, 2018) sited at the University of Oxford (Jolley and Maiden 2010).

Table 1

Strains Sequenced in This Study, Origin of Specimen and Genome Sequencing Data

Strain	Year of Isolation	Sequence Type (ST)	Origin of Specimen	Genome Size (kb)	Contigs (#)	Coverage (fold)	N50 (kb)	N90 (kb)
HIAE_PA01	2011	455	Blood	6,764	77	68	266	68
HIAE_PA02	2012	381	Blood	6,459	44	175	322	97
HIAE_PA05	2013	498	Blood	6,803	61	163	245	92
HIAE_PA06	2014	231	Blood	6,328	44	92	622	98
HIAE_PA07	2015	252	Blood	6,970	124	177	212	41
HIAE_PA08	2015	new	Bronchoalveolar lavage	6,443	30	175	502	142
HIAE_PA09	2015	245	Blood	6,575	58	87	397	94
HIAE_PA10	2015	2235	Bronchoalveolar lavage	6,527	106	89	227	62
HIAE_PA11	2015	244	Bronchoalveolar lavage	6,807	51	181	328	122
HIAE_PA12	2015	10	Blood	6,482	40	83	305	107
HIAE_PA13	2015	244	Blood	6,956	113	83	234	59
HIAE_PA14	2015	1290	Bronchoalveolar lavage	6,320	23	172	1,331	210
HIAE_PA15	2016	1993	Blood	6,304	37	189	521	99
HIAE_PA16	1994	606	Bronchoalveolar lavage	6,869	80	89	306	61
HIAE_PA17	1995	235	Tracheal secretion	6,735	59	114	273	94
HIAE_PA18	1995	348	Bronchoalveolar lavage	6,765	81	213	220	77
HIAE_PA19	1995	235	Tracheal secretion	6,847	88	186	228	83
HIAE_PA20	1997	235	Tracheal secretion	6,812	79	258	237	85
HIAE_PA21	1998	253	Tracheal secretion	6,862	97	161	237	68
HIAE_PA22	1995	190	Surgical secretion	6,500	50	184	316	98

Besides the 20 genomes sequenced in this study, we used all closed genomes of P. aeruginosa strains that were available from GenBank in December 2017 (in total, 124 available genomes), and 79 genomes of ST235 strains available in multiple contigs for phylogenomic analysis. Gene prediction and annotation of all genomes were done with RAST (Aziz et al. 2008). For the phylogenomic analysis, the core genome was aligned using Parsnp, a program that is part of the Harvest software package (Treangen et al. 2014). All sequences near the end of contigs (gaps) were excluded for this analyses in order to avoid using false SNPs for phylogenetic reconstruction (default settings in Parsnp). Phylogenetic trees were built in Mega v7 and visualized with iTOL (Kumar et al. 2016; Letunic and Bork 2016). To determine orthologous genes among the 20 sequenced P. aeruginosa strains we used Proteinortho (Lechner et al. 2011). As parameters the default settings of the program were used: BlastP e-value cutoff, 1e−05; minimum identity of best blast hits, 25%; minimum coverage of best blast hits, 50%. Genome comparison figures were generated using BLAST Ring Image Generator (BRIG) and based on BlastN analyses (Alikhan et al. 2011). Proteinortho and BRIG were used in combination to determine the accessory genome. For functional predictions and pathway analyses we used KEGG mapper (http://www.genome.jp/kegg/mapper.html; last accessed June 19, 2018).

Results and Discussion

WGS of 20 epidemiologically unrelated strains collected over a 20 year-period resulted in final genome assemblies with 23 to 124 contigs (average: 67) per strain; N50 and N90 values ranged from 212 to 1,331 kb (average: 366 kb) and 41 to 210 kb (average 93 kb), respectively (table 1). The GC-content of the strains varied between 65.7% and 66.5%, and the genome size varied between 6,304 kbp and 6,970 kbp. The presence of different STs among the 20 strains illustrates the heterogeneity of the cohort.

The pan-proteome (all orthologous, core genome-based CDS and all CDS of the accessory genome) of the 20 strains compared with reference strains PAO1 and UCBPP-PA14 comprises 7,874 CDS. The 22 strains share 4,844 orthologous CDS, representing the core functions as well as a high array of two-component systems (and proteins involved in quorum sensing, as predicted by KEGG pathway analyses (data not shown). In total, 117 CDS are strain-specific and 626 proteins are only found in maximal two genomes, representing the flexible genome (supplementary table S1, Supplementary Material online). Besides phage-related proteins, mobile elements, prophages and restriction-modification systems, the flexible genome encodes transport and resistance-related functions as well as proteins for the degradation of toxic substances (biphenyl dioxygenase, lactoylglutathione lyase, acriflavine degradation). The three ST235 strains harbored the most diverse set of acquired resistance-associated genes, including the insertion sequences int1, tn7, and tn21 and also contained aac(6’)Ib-cr encoding an aminoglycoside modifying N-acetyltransferase reported to confer resistance to fluoroquinolone (supplementary table S2, Supplementary Material online). The gene aac(6’)Ib-cr has recently been reported for the first time in nosocomial P. aeruginosa strains isolated in Brazil (Araujo et al. 2016). Our results indicate, however, that the gene was already present in P. aeruginosa ST235 strains collected as early as 1995 in the country. This plasmid-mediated quinolone resistance determinant was described in P. aeruginosa ST235 for the first time in 2008 (Libisch et al. 2008). Interestingly, the ST235 strains carrying this gene in this study also carried intI suggesting that a class 1 integron may be harboring this resistance gene as in ST235 strains from Hungary and Serbia (Libisch et al. 2008). This finding highlights how WGS of clinical isolates may contribute to the understanding of resistance markers spread, worldwide, through high-risk strains.

Noteworthy was also the widespread presence of the acquired blaOXA-50 and blaPAO genes, contributing to β-lactam resistance, and of the gene aph(3’)IIb, conferring resistance to aminoglycosides, in strains from both studied periods. The presence of genes encoding effectors secreted via the type III secretion system (exoS, exoT, exoU, and exoY) in clinical isolates of P. aeruginosa is variable. Here the three ST235 isolates were exoU-positive; three other STs, including a not previously reported one, also carried exoU. This gene encodes a major virulence factor, a phospholipase, which has been associated with poor clinical outcomes (Finck-Barbançon et al. 1997; Sato et al. 2003). Besides gene acquisition, the selection of chromosomal mutations is also a frequent occurrence associated with multidrug resistance in P. aeruginosa. Resistance to fluoroquinolones is mostly associated with the Quinolone-Resistance-Determining Region (QRDR) (Jaillard et al. 2017; del Barrio-Tofino et al. 2017). The three ST235 strains showed acquired mutations in QRDR, as also recently described in other Brazilian isolates (Araujo et al. 2016). The identification of antibiotic resistance markers in the genome correlates and often broadens observations obtained by traditional phenotypic assays (Jaillard et al. 2017). Taken together, WGS helped in the identification of persistent resistance and virulence markers in P. aeruginosa strains within a clinical setting that may reflect antibiotic policies.

As previously reported, P. aeruginosa population consists of three clades. Since most sequenced strains to date are isolates from human infections there could be a bias towards clinical isolates, but there is currently no convincing explanation for this specific population structure of P. aeruginosa (Hilker et al., 2015; Wiehlmann et al. 2015). It has been suggested that different environment habitats/niches and/or different host associations led to this evolutionary bifurcation (Freschi et al. 2015; Hilker et al. 2015; Wiehlmann et al. 2015). A phylogenetic analysis based on core genome-based SNPs separated the 20 P. aeruginosa strains sequenced in this study into these previously reported groups (fig. 1). One strain, HIAE_PA16, is a hybrid strain with genomic regions from both clades; it belongs to the previously described group 3 that contains only few strains sequenced to date (Freschi et al. 2015). Additionally, this analysis showed that most of the older strains, isolated between 1994 and 1998, are phylogenetically distinct from strains isolated more recently (between 2011 and 2016). We were not able to find other reports that indicate this group shift, but we are also not aware of a study that has investigated nosocomial P. aeruginosa strains collected over a larger time period at the phylogenomic level.

. 1.

—(A) Phylogenomic comparison of the 20 P. aeruginosa strains isolated in a Brazilian hospital over a time period of two decades. The phylogenetic relation was reconstructed by a core genome alignment and the comparison of high-quality SNPs. Two reference strains with closed genomes were added, PAO1 and UCBPP-PA14. The three ST235 strains (HIAE_PA17, HIAE_PA19 and HIAE_PA20) are closely related, as well as the two ST244 strains (HIAE_PA11 and HIAE_PA13), and the two ST253 strains (HIAE_PA21 and the reference strain UCBPP-PA14). There is also some relatedness between the strains HIAE_PA05 (ST498) and HIAE_PA07 (ST252). The branch length indicates the individuality of each strain. (B) Phylogenomic analysis of all previously sequenced strains of P. aeruginosa and the 20 strains isolated from Brazil. Only strains with complete genomes were analyzed, altogether 124 strains, together with the here sequenced 20 strains (in red, strains isolated between 1994 and 1998; in blue, strains isolated between 2011 and 2016). The analysis revealed that most strains, except three (belonging to group 3, nodes in black), grouped into two clades: 45 strains formed group 2 (clade in orange), and 96 strains formed group 1 (clade in purple). (C) Core genome SNP analysis of ST235 strains and the phylogenomic reconstruction. Phylogenetic analyses based on SNPs located in the core genome of all publicly available ST235 strains. In red, the three ST235 strains sequenced in this study; in orange, strains from the same clade (they are all strains isolated in South America). The respective clade is highlighted in green. (D) Visualization of genome comparison of the 20 Brazilian strains. A BRIG genome comparison analysis was performed with all strains isolated in this study and the two strains PAO1 and UCBPP-PA14. The analysis was carried out with HIAE_PA17 as the reference genome. Differences regarding the flexible genome can be seen.

In order to further explore this clade distribution of the here sequenced P. aeruginosa strains, we analyzed all available P. aeruginosa strains for which a closed genome sequence was available (fig. 1). The total number of SNPs in the determined core genome (62% of the PAO1 reference genome) was 185,860. The analysis showed that almost all P. aeruginosa strains (98%) clustered into the above-mentioned group 1 and group 2, and three strains were considered possible hybrids of group 1 and group 2 genomes, belonging to the previously described group 3 (Mathee et al. 2008). There were few clade-specific genes present (data not shown). The presence of the genes exoS and exoU genes has been reported as mutually exclusive (Feltman et al. 2001). In fact, they were found only in group 1 and group 2 genomes, respectively. The gene exoS encodes the toxin exoenzyme S, which is an actin-ADP-ribosylating toxin, delivered by the type III secretion system (Engel and Balachandran, 2009).

We decided to have a closer look on group 2 since it contained strains belonging to ST235, a predominant global, high-risk clinical clonal lineage (Treepong et al. 2018). We wanted to know how closely related the three here sequenced ST235 strains are and how they relate to previously identified ST235 strains. In total, 79 genomes, for which metadata were available were used (Treepong et al. 2018). The respective strains originated from various geographic locations and were isolated between 1991 and 2014. The ST235 strains shared a core genome of 72% with a total number of 4,302 SNPs. The three Brazilian strains sequenced here clustered in close vicinity, revealing that they are closely related to each other and the closest relatives of these three strains were strains either isolated in Brazil or Argentina (fig. 1). In order to visualize the core and the flexible genome of ST235 strains, a BRIG analysis was performed, using HIAE_PA17 as the reference genome (fig. 1). Strain HIAE_PA17 contained an 85 kb genomic island (GI), inserted next to a tRNA (tRNA-Gly) (fig. 1). This GI has signatures of an integrating conjugative element (ICE), similar, but not identical, to other P. aeruginosa genomic islands (PAGIs), such as PAGI-2 and -3 (Kung et al. 2010); it encodes among many other functions a type IV secretion system (GI-type T4SS). There are two other GIs in this strain, carrying a GI-type T4SS; one of which is inserted next to a tRNA (tRNA-Lys) and contains the exoU gene. The latter is also present in the other two ST235 genomes. Strain HIAE_PA20 possesses a strain-specific 36 kb prophage, and a second 37 kb prophage that is also absent in the other sequenced strains, except for HIAE_PA19.

Taken together, these results demonstrate the closer relationship between ST235 strains isolated in Brazil and Argentina when compared with ST235 strains isolated in other geographic locations, a noteworthy report since there is no epidemiological link between these isolates. However, the SNP differences in the core genome as well as the different genomic islands in the flexible genome underline the individuality of the three ST235 strains isolated in the same clinical setting and possibly emphasize the capacity of ST235 strains to quickly adapt and modify its genome according to its environment. In fact, WGS is becoming more accessible and the analysis of entire genomes may provide enough resolution for the discrimination of highly related lineages such as ST235 in routine screenings.

In conclusion, analysis of our strain cohort corroborates data of recent studies that have presented two major clades for the classification of P. aeruginosa population. Reasons for this distribution remain to be found. Each strain showed a distinct evolutionary history as judged from the many strain-specific core genome-based SNPs and the presence of strain-specific GIs. This illustrates the high diversification of the species P. aeruginosa in this clinical setting, including the three ST235 isolates. The three ST235 strains were highly similar to other ST235 isolates from Brazil and Argentina, but not with isolates from other parts of the world, a new insight into epidemiologically unrelated ST235 strains. WGS also allowed for the identification of resistance markers that have persisted in P. aeruginosa genomes during the past two decades, despite dissimilar antibiotic policies during this period.

Supplementary Material

Supplementary data are available at Genome Biology and Evolution online.

Click here for additional data file.