Evolution and host-specific adaptation of Pseudomonas aeruginosa

Abstract

The major human bacterial pathogen Pseudomonas aeruginosa causes multidrug-resistant infections in people with underlying immunodeficiencies or structural lung diseases such as cystic fibrosis (CF). We show that a few environmental isolates, driven by horizontal gene acquisition, have become dominant epidemic clones that have sequentially emerged and spread through global transmission networks over the past 200 years. These clones demonstrate varying intrinsic propensities for infecting CF or non-CF individuals (linked to specific transcriptional changes enabling survival within macrophages); have undergone multiple rounds of convergent, host-specific adaptation; and have eventually lost their ability to transmit between different patient groups. Our findings thus explain the pathogenic evolution of P. aeruginosa and highlight the importance of global surveillance and cross-infection prevention in averting the emergence of future epidemic clones.

Figure 1. The emergence of epidemic clones of Pseudomonas aeruginosa.

(A) Geographical location of the whole genome sequenced P. aeruginosa isolates obtained from patients, animals, and environment analysed in this study (n = 9,573). Number of samples from each location indicated by the size of blue dot. (B) Cumulative number of isolates across P. aeruginosa clones (defined by clustering genomes using the unweighted pair group method with arithmetic means; see Supplementary Methods), arranged by ascending number of genomes per clone and stratified into epidemic (n ≥30 isolates/clone; red), non-unique (1 < n < 30 isolates/clone; light brown), and unique (n = 1 isolate/clone; blue) groups. (C) Left: Maximum likelihood phylogenetic tree generated from genomes of all study isolates (major epidemic clones labelled in red). Right Bar plot representing the number of cities where each epidemic clone was found, coloured by continent. (D) Estimated date of first population expansion of 21 epidemic clones (predicted by Bayesian inference using BEAST (29)) with graph showing median and interquartile range (IQR; boxplots), 1.5 times IQR range (whiskers), and data points outside this range (black points). (E) Pangenome graph analysis of ancestral representatives of epidemic clones (n = 21) and sporadic clones (n = 80), constructed using Panaroo (39), where nodes represent clusters of orthologous genes and two nodes are connected by an edge if they are adjacent on a contig in any sample from the population, define gene gain events associated with the emergence of epidemic clones (described in detail in Figure S5) with genes highlighted that are involved in transcription (blue), defense mechanisms (purple), and inorganic ion transport and metabolism (yellow).For illustration purposes, the graph has been ordered against the genome of P. aeruginosa PAO1. Inset: magnified section of the pangenome graph is shown to illustrate node and edge structure.

Figure 2. Variable intrinsic host preference of epidemic P. aeruginosa clones.

(A) Proportion of infections caused by epidemic clones (labelled by their majority multi-locus sequence type, ST) in cystic fibrosis (CF; red) and non-CF (blue) patients. (B) UMAP projection of transcriptomes from representative isolates of epidemic clones (25), colour-coded by the CF affinity of each clone. Expression data were pseudo-aligned to strain-specific gene indices to produce estimates of gene transcript abundance. (C) Transcriptome-wide association of gene expression with CF affinity. Transcript abundances were modelled as a response to the proportion of CF infections caused by each epidemic clone using a negative binomial generalised linear model. Volcano plot visualization of the Log2-fold expression change with CF proportion for every gene in the 99% core genome of Pseudomonas aeruginosa (center). Genes with an adjusted p-value of less than 0.05 and a log2 fold change less than -0.5 were coloured in green, genes with a log2 fold change greater than 0.5 were coloured in red. The coefficients for gene models were assessed using the Wald test (FDR = 0.05). Normalized expression counts vs CF proportions per epidemic clone with a trendline for the two genes with the lowest and highest log2 fold change, respectively, are shown above (top left/top right). Bulk RNA seq data was analysed from 241 clinical isolates of epidemic clones (25) included in our strain collection. (D) Survival of epidemic clones within wildtype (WT) or isogenic F508del knock-in THP1 macrophages at 2 and 4h post infection, expressed as fold change from 1 hour post infection showing median and interquartile range (IQR; boxplots), 1.5 times IQR range (whiskers. Experiments (carried out at least in duplicate) were performed by exposing THP1 macrophages to pooled isolates of 51 clinical isolates at a multiplicity of infection (MOI) of less than 1. Viable bacteria were isolated from macrophages at time points indicated and grown on solid media. Isolate abundance was quantified using sequence-based deconvolution. Strains with less than 1% abundance at the 1h time point were excluded from the analysis. A difference in the abundance of ST27 strains vs ST111 and ST235 strains at the 4h timepoint was assessed using a two-tailed t-test. * p-value < 0.05, ** p-value < 0.01.

Figure 3. Activation of the DksA1 regulon contributes to Cystic Fibrosis host preference of P. aeruginosa clones.

(A) Volcano plot visualisation of the Log2-fold expression change with CF proportion for genes positively controlled (red) and negatively controlled (green) within the DskA1 regulon as defined by Fortuna et al. (44). Bulk RNA seq data was analysed from 241 clinical isolates of epidemic clones (24) included in our strain collection. (B) DksA1 promotes survival of P. aeruginosa within CF macrophages. Viable intracellular P. aeruginosa (quantified through enumeration of cell-associated colony forming units; CFU) were measured at 1h and 4h post infection of differentiated wildtype (WT) and isogenic F508del homozygous knockin (CF-F508del) THP1 cells with wildtype (blue), isogenic DskA1-DskA2 double knockout (ΔDksA1,2; pink), and knockout complemented with DksA1 ((ΔDksA1,2::DksA1; yellow) P. aeruginosa PAO1. Data (mean ± SEM) are representative of at least three independent experiments performed in at least triplicate. *** p < 0.001; ns not significant (two-tailed Student’s t-test). (C) (B) Top: Cartoon of zebrafish (created with BioRender.com) illustrating injection site for GFP-labelled fluorescent P. aeruginosa. Bottom: Representative fluorescence and DIC images of whole infected zebrafish larvae at 1 day post-infection (Scale bar: 150 μm; the labelled yolk sac is autofluorescent). (D) Survival analysis of control (top) and cftr morphant (cftr MO; bottom) zebrafish larvae infected intravenously (250-350 CFU) with P. aeruginosa PAO1 wildtype (blue), ΔDksA1,2 knockout (pink), and ΔDksA1,2::DksA1 complemented (yellow) fluorescent strains plotted as the percentage of surviving animals over 6 days (average of 2 independent experiments; n = 30-38 fish for each condition); *** p < 0.001 (Mantel-Cox Log-rank test). (E) Viable P. aeruginosa in zebrafish larvae at Day 1 post infection with P. aeruginosa PAO1 wildtype (blue), ΔDksA1,2 knockout (pink), and ΔDksA1,2::DksA1 complemented (yellow) fluorescent strains (plotted as mean ± IQR colony forming units (CFU) per fish of at least 3 independent experiments; n = 15-20 larvae per condition. *** p < 0.001; ns not significant (two-way ANOVA with Tukey’s post-test). (F,G) Control and cftr morphant zebrafish larvae with mCherry-labelled macrophages (Tg(mpeg1:mcherry-F)ump2 (45)) were intramuscularly infected with 250-350 GFP-labelled P. aeruginosa PAO1 wildtype, ΔDksA1,2 or ΔDksA1,2::DksA1 strains) and the infection tracked using real-time intravital confocal microscopy. (F) Representative 3D reconstruction of confocal imaging showing macrophages (red) and automatic classification of extracellular (grey) and intracellular (green) P. aeruginosa (Scale bar 10 μm). (G) Quantification of the number of infected macrophages at the site of injection (left) and the level of intracellular bacterial load (calculated by the volume of bacteria-associated fluorescence observed within each macrophage) at 6 hours post infection with P. aeruginosa PAO1 wildtype (blue), ΔDksA1,2 knockout (pink), and ΔDksA1,2::DksA1 complemented (yellow) fluorescent strains. Mean ± IQR of at least 54 cells per condition (from n = 4-6 larvae) recorded from 2 independent experiments. ** p < 0.01; *** p < 0.001; ns not significant (two-way ANOVA with Tukey’s post-test).

Figure 4. Host-specific pathoadaptation of P. aeruginosa

(A) Manhattan plot showing nominal p values (plotted as -Log10) from genome-wide mutational burden test across all genes in P. aeruginosa PAO1. Significance was assessed using a Poisson test comparing the expected and observed number of mutations in each gene accounting for the proportion of genomes that gene was found in the pan-genome (FDR = 0.1; genes with a significant mutational burden, termed pathoadaptive, shown in black, others in grey). (B) UMAP projections of host adaptation of isolates (based on acquired mutations in pathoadaptive genes) colour-coded by (left) number of pathoadaptive mutations and (right) type of infection (centroids denoted by larger dots). Isolates without any pathoadaptive mutations were removed from the analysis. (C) Protein-protein interaction network for the pathoadaptive genes extracted from the STRING database (only main connected component shown, full graphs shown in Fig. S12; (56)). Genes are shown as nodes which are connected by an edge if they had an interaction reported in STRING (confidence > 0.7). Top: To estimate host-specific pathoadaption, the number of cystic fibrosis (CF) vs non-CF mutations (determined by stratifying mutations in pathoadaptive genes on terminal branches by the infection type of isolates) were compared using a Fisher exact test (FDR = 0.1) and expressed as an odds ratio for each gene. Bottom: Gene nodes were colour-coded by class of functional annotation (based on overrepresented pathways using Gene Ontology (89) biological process enrichment analysis with TopGO (57) among CF: transmembrane transport and fatty acid biosynthesis, and non-CF: transcriptional regulation and chemotaxis).

Figure 5. Evolutionary trajectories of P. aeruginosa during pathoadaptation.

(A) Normalised frequency of mutations over evolutionary time in specific pathoadaptive genes. The trajectories of the 50 most commonly mutated genes were manually assigned to one of 5 classes (Figure S15), based on the shape of their mutation frequency curves (relative size of each class and representative examples (with trendlines from locally-weighted smoothing) shown). (B) The relative transmissibility and host-specific adaptation of pathoadaptive genes was calculated. To estimate host-specific pathoadaptation, the number of cystic fibrosis (CF) vs non-CF mutations (determined by stratifying mutations in pathoadaptive genes on terminal branches by the infection type of isolates) were compared using a Fisher exact test (FDR = 0.1) and expressed as an odds ratio. To assess the transmissibility of pathoadaptive changes, the number of mutations that had been observed in at least two isolates were compared with mutations that had only been observed once using a Fisher exact test (FDR = 0.1). Genes were colour-coded if showing significant host-specific adaptation (blue), changes in transmissibility (purple), or both (pink). Genes with zero or infinite odds ratio not shown. (C) Functional annotation of pathoadaptive genes associated in (top) host-specific adaptation and (bottom) changes in transmissibility. (D) The number (top) and proportion (bottom) of transmission links across a range of pairwise SNP thresholds occurring between CF to CF (red), CF to non-CF (yellow), and non-CF to non-CF (blue) individuals (data were down-sized to contain equal numbers of CF and non-CF infections). (E) Transmission clusters involving patients with CF (red), non-CF (blue), or unknown status (white). Nodes representing isolates were connected by edges if pairwise SNP distances were 26 SNPs or less. This cut-off represents the 95th percentile of the within-host genetic diversity analysed in 81 patients.