P. aeruginosa 6378/09 strain (NMI collection number), hereinafter referred to as clinical strain and P. aeruginosa ATCC 27853 control (reference) strain, obtained from American Type Culture Collection, were used for proteomic analysis. It is wildtype, imipenem-sensitive strain, hereinafter referred to as control strain. The P. aeruginosa clinical strain was isolated from an oncology patient from an abdominal swab. It is an ESBL-negative, XDR strain, resistant to all important in therapy antibiotics (amikacin, aztreonam, cefepime, ceftazidime, ceftazidime/avibactam, ceftolozane/tazobactam, ciprofloxacin, colistin, imipenem, imipenem/relebactam, levofloxacin, meropenem, meropenem/vaborbactam, piperacillin/tazobactam, tobramycin) except colistin (MIC=2 mg/L). Resistance to imipenem is associated with VIM-2 variant of metallo-beta-lactamase (MBL), and the strain belonged to world epidemic clone ST235, as previously described (Urbanowicz et al., 2021). ST235 has been the most frequent MDR/XDR P. aeruginosa clone globally, associated with VIM/IMP-like MBLs (Oliver et al., 2015) and VIM-2 are the most common variant of MBL globally (Mack et al., 2025).

The genomes of P. aeruginosa ATCC 27853 (Cao et al., 2017; BioProjectPRJNA377172, Accession number: CP011857) and NMI clinical strain ST235 (Urbanowicz et al., 2021; BioProject PRJNA701400, BioSample: SAMN17864755) were recruited from Genbank and compared using Bakta (ref: https://www.microbiologyresearch.org/content/journal/mgen/10.1099/mgen.0.000685), steps followed are available at https://github.com/animesh/pseudomonas_align_strains?tab=readme-ov-file#exact-commands-i-ran-for-bakta-in-this-project.

Minimum inhibitory concentration (MIC) of imipenem (Merck KGaA, Darmstadt, Germany) was determined for both strains by the broth microdilution method according to standard ISO 20776-1 (ISO 20776-1, 2020). MIC (mg/L) ranges used in the tests were 1 – 512. The MIC was defined to be the lowest concentration at which no visible growth (complete inhibition of growth) of bacteria could be observed after incubation for 18 +/- 2 hours as recommended by the European Committee on Antimicrobial Susceptibility Testing (European Committee on Antimicrobial Susceptibility Testing, 2024).

The bacteria were grown in cation-adjusted Mueller-Hinton broth (Becton Dickinson, Sparks, MD, USA) at 37°C with shaking (120 rpm) to optical density OD₆₀₀ = 1.0. For antibiotic challenge, imipenem concentration of 0.1x MIC, 0.25x MIC, 0.5x MIC and 1x MIC of the particular strain was used in a volume of 5 ml. Bacterial cultures without the antibiotic were considered as untreated controls (see Figure 1A ). Challenged bacterial cultures were then incubated with shaking (120 rpm) for 5 hours at 35°C +/- 2°C. Each experiment was performed three times for both clinical and control strains (3 biological replicates) with two technical replicates for each condition (antibiotic concentration). After incubation the bacterial cell pellets were harvested in pre-cooled Falcon tubes by centrifugation (5 min at 3,000 × g below 0°C). Then, the cell pellets were washed twice with ice-cold 1 x PBS (Sigma-Aldrich, Life Science, St. Louis, MO, USA). Finally, cell pellets were collected by centrifugation (5 min at 3,000 × g below 0°C), snap-frozen, and stored at −80°C until further testing (Liu et al., 2016; Hashemi et al., 2019; Jongers et al., 2021).

Proteins were quantified by processing the timsTOF-pro DDA MS-data using MaxQuant v2.0.3.0 (Tyanova et al., 2016). Namely, the following search parameters were used: enzyme specified as trypsin with a maximum of two missed cleavages allowed; minimum peptide length of 7; acetylation of protein N-terminal, oxidation of methionine, and deamidation of asparagine/glutamine as dynamic post-translational modification. These were imported in MaxQuant which uses m/z and retention time (RT) values to align each run against each other sample with a minute window match-between-run function and 20 mins overall sliding window using a clustering-based technique. These were further queried against the reference-proteome of P. aeruginosa downloaded from UniProt (UniProt Consortium, 2022) in May 2022 and MaxQuant’s internal contaminants database using Andromeda built into MaxQuant. Both Protein and peptide identifications false discovery rate (FDR) was set to 1%, only unique peptides with high confidence were used for final protein group identification. Peak abundances were extracted by integrating the area under the peak curve. Each protein group abundance was normalized by the total abundance of all identified peptides for each run and protein by calculated median summing all unique peptide-ion abundances for each protein using label-free quantification (LFQ) algorithm (Cox et al., 2014) with minimum unique peptide(s) ≥ 1. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE (Perez-Riverol et al., 2022) partner repository with the dataset identifier PXD055744 (Username: reviewer_pxd055744@ebi.ac.uk Password: XNtCwVbI4AIs).

Having obtained the label-free quantification (LFQ) intensity values for both the P. aeruginosa clinical strain and control strain, across the five different concentrations (as depicted in Figure 1A ), we performed a series of data filtering, processing, and aggregation steps, including removal of reverse hits and potential contaminants to ensure the integrity of our data. In the case of multiple UniProt protein IDs per hits, the one with the highest curation status were used in the downstream analysis when relevant. This led to a total of N=3811 proteins in our analysis, accounting for approximately 60% of the chromosomal genes in an average P. aeruginosa strain (Mosquera-Rendón et al., 2016).

For each strain and antibiotic concentration, we had three biological replicates ( Figure 1A ). We performed a standard t-test using the log2-transformed LFQ (log2LFQ) values for each protein ID to obtain the associated p-values. These p-values were then corrected using the Benjamini-Hochberg False Discovery Rate (FDR) procedure. This was done for all experiments with non-zero MIC, using the 0.0x MIC as a control for each strain. This means that we compared LFQ values and performed t-tests within strains only. Alongside the t-test, p-value and the FDR values, we used the median value for LFQ intensities where applicable.

When comparing LFQ intensity values between experiments, we primarily used fold change relative to the control probe (no antibiotic) to describe the difference in log2 space, expressed via log2LFQ values. Proteins with LFQ intensity values of zero for both the control probe (no antibiotic) and 1.0x MIC experiments were excluded. Proteins with values of zero in one experiment and non-zero in another were assigned a large, finite fold-change value. When comparing across strains, we used the median log2LFQ fold changes between the strains, as these values were already normalized with respect to the 0.0x MIC control experiment for the respective strains.

Unless otherwise specified, the intra-strain protein LFQ intensity fold-change cutoff was set at 1.5, or approximately +/- 0.5849 after log2 transformation. In the individual protein volcano plots an uncorrected p-value cutoff of 0.1 was employed, with values above colored grey and values below colored light red or blue, while FDR values below 0.1 are furthermore indicated by dark red or blue color. These thresholds are not particularly conservative, which is deemed appropriate given the explorative nature of the study. The individual protein volcano plot x-axis and y-axis ranges were fixed to improve visual comparability, and log2 fold changes were visually capped at +/- 5 alongside p-values below 1e-5.

To visualize the overall features of the data set, we calculated and used a Principal Component Analysis (PCA), using standard scaling, and a correlation matrix (heatmap) of Spearman’s correlation coefficients between experiments. This provided an overview of the structure and relationships in the data set ( Figures 1B, C ).

To discern significant systemic alterations in the proteomes within and across strains and antibiotic concentrations, we conducted a Gene Ontology (GO) term enrichment analysis. This method is set-based, meaning that groups of proteins are compared rather than individual proteins. Statistical tests are performed by comparing a subset to the full background set. The GO term enrichment analysis aims at finding enriched functional, process-related, and compartmental annotations associated with each protein in the subset. Collectively, numerous individually insignificant protein-level changes can be linked to statistically significant features on aggregate, cellular level.

The GO term enrichment analysis is agnostic to the subset selection method. Therefore, we employ both FDR-value, uncorrected p-values, and no p-value thresholds when selecting proteins, in conjunction with fold-change thresholds, at various parts of the analysis. We performed statistical tests based on hypergeometric distributions on the subset selection, and subsequently corrected these using the Benjamini-Hochberg FDR procedure, setting the FDR-value threshold to 0.1. Both over-representation and under-representation were accounted for in the analyses. The choice of background, i.e., the full set against which the subset is compared, is crucial for obtaining accurate results. For this, we used all mapped and processed proteins with well-defined t-tests based on the three biological replicates.

We collected GO term information for P. aeruginosa using UniProt, release: 2024_03 (UniProt Consortium, 2023), and mapped the GO term IDs onto the Gene Ontology Knowledgebase using AmiGO 2, version: 2.5.17 (Ashburner et al., 2000), aligning the information with the UniProt IDs. Out of 3811 proteins, 3253 were associated with GO terms, and the remaining were assigned a placeholder value to ensure their inclusion in the analyses, since the non-annotated proteins may be of particular interest for further study.

In the results, we visualized the enriched GO terms using volcano plots. The x-coordinate of each GO term represents the log2 value of its enrichment (blue indicates under-representation and red indicates over-representation), and the range is fixed and visually capped at +/- 2.5. The y-coordinate represents the FDR values associated with each GO term (not to be confused with the protein-level FDR-values). The area of the GO-term circles is proportional to the number of proteins associated with the GO term in the subset. All terms with a GO-term FDR value > 0.1 are indicated in grey and are excluded from further analysis.

Data analyses and visualizations were performed using Python 3.11, along with standard modules such as numpy, pandas, scipy, matplotlib, and seaborn, among others. The t-tests and dose-correlation statistics were conducted using R.

Comparative genomes analysis revealed that both strains demonstrate a significant degree of similarity, with nearly 4,000 potential gene products shared between them. The fundamental functions of Pseudomonas aeruginosa remain preserved. Notwithstanding the phylogenetic divergence, there exists considerable overlap alongside a common pool of potential target proteins. However, the control strain lacks the distinctive set of resistance factors present in ST235, which confer a strategic advantage in clinical settings ( Supplementary Figure S2 ) along with 573 and 835 hypothetical protein in ATCC 27853 and ST235 respectively ( Supplementary Table T2 ).

Exploration of the overall trends in the data, as visualized in Figures 1B, C , suggests that the protein abundance in the clinical strain remains stable under the influence of imipenem, unlike the control strain. The PCA plot ( Figure 1B ) reveals that the largest variance between the experiments is attributed to the strains themselves (PC1, x-axis). The second largest factor of variance (PC2, y-axis) corresponds to the increasing imipenem concentration, as evidenced by the gradual transition from control probe (no antibiotic) in the top-right to 1.0x MIC in the bottom-right for the control strain. The heatmap ( Figure 1C ) displays a highly intercorrelated block in the top-left corner, representing the clinical strain experiments. The bottom-right corner shows a slight but clear gradient towards covariance among experiments with similar imipenem concentrations. Notably, there is high similarity between the strains, with correlation coefficients around 0.85 for the control probes (A1 and B1) and 0.80 at 1.0x MIC (A5 and B5).

Based on the overall results from Figures 1B, C , we expect the significant fold changes for individual proteins across MIC concentrations to be more pronounced for the control strain than the clinical strain. The volcano plots in Figure 2 confirms this prospect. Here, we have plotted the N=3811 proteins with a fold change threshold of 1.5 for both the control ( Figure 2A ) and clinical ( Figure 2B ) strains at 1.0x MIC, as compared with a t-test to 0.0x MIC (see Materials and Methods for details).

In these plots, an increase in fold change from 0.0x MIC to 1.0x MIC, i.e. up-regulation, is depicted by red dots, while a decrease in fold change, i.e. down-regulation, is represented by blue dots. Each dot corresponds to a protein. For both strains, we display the proteins with p-values above 0.1, and fold change below the threshold, as grey, while light red and blue denote p-values below 0.1 and full (darker) red and blue depict the proteins with FDR corrected p-values below 0.1. As expected, the control strain shows a range of significant FDR values, while the clinical strain only has 8 proteins fulfilling the criteria: AmpC (P24735), PA0868 (Q9I573), DDAH (Q9I4E3), Tse4 (Q9I069), PA2547 (Q9I0T7), PA0207 (Q9I6T1), VreR (G3XCU5) and PhzB2 (Q9S508). These proteins (depicted in the top corners of Figure 2B ) are further explored in Figure 3 and Supplementary Table S1 and detailed results will be presented in the following subsection.

For the normalized 1.0x MIC control strain experiments, displayed in Figure 2A , there are n=601 proteins with FDR value below 0.1 which show a fold change greater than 1.5 (depicted in dark red and blue), and n=1219 proteins adhere to the fold change threshold alone. The respective numbers for the clinical strain are n=8 (where n=79 have an uncorrected p-value below 0.1) and n=358, as visualized in Figure 2B .

As mentioned in the previous section, 8 proteins passed the FDR value threshold of 0.1 for the clinical strain, making them resistome candidates which are significantly involved in the antimicrobial resistance in P. aeruginosa. These proteins are shown in Figure 3 , where Figures 3A-D show an overall down-regulation of the specific protein in the clinical strain (in green) for increasing MIC, Figures 3E-F show an overall up-regulation for increasing MIC, while Figures 3G-H show a particularly large difference between the strains, i.e., where the LFQ intensity in the clinical strain is much lower than the control strain. The control strain (in orange) shows an overall trend towards up-regulation for these 8 proteins, aside from the proteins in Figures 3C, D . The standard deviations, indicated by shaded green for the clinical strain and shaded orange for the control, vary. Large standard deviations, such as for PA0868 (Q9I573) in Figure 3D , indicate high variability in protein expression, while the high-LFQ intensity protein AmpC (P24735) in Figure 3F displays a clear signal, suggesting consistent protein expression across replicates. The dose-correlation statistics were calculated for all proteins for both strains, but since the number of measurements per protein is low, the associated p-values are poor, potentially making the correlation values misleading. Consequently, the specific dose correlation values will not be discussed here, but please refer to Supplementary Tables S3 and S4 for details on clinical and control strains respectively.

The proteins in Figure 3 are described in more detail in Supplementary Table S1 , where a description and associated gene ontologies are added (see Materials and Methods for details). These findings and their implications for understanding antimicrobial resistance in P. aeruginosa will be discussed further in the Discussion section.

To better understand the biological implications of the significant up- and downregulated protein subsets (highlighted in red and blue in Figure 2 ), we conducted a GO term enrichment analysis as outlined in the Materials and Methods section. This analysis allows us to assess the biological relevance of our protein subsets, indicated by the enrichment of biological traits, when compared to the original full set of N=3811 proteins. By comparing high-fold-change proteins between 1.0x MIC and control probe (no imipenem), we can identify the biological functions and processes that the organism upregulates or downregulates. This comparison also provides p-values, which can furthermore be FDR corrected, indicating the statistical significance of these enrichments. I.e., there are two levels of (corrected) p-values at this stage: (i) protein-level p-values, based on t-tests for individual proteins, as depicted in Figure 2 , and (ii) protein-subset level p-values, where enrichments of protein annotations in the subset compared to the full, original set are assessed via hypergeometric processes. Thus, a random selection of proteins, which does not reflect real biology, is expected to meet the protein-subset level, GO-term FDR threshold of 0.1, thus reducing the overall likelihood of false positive GO-term enrichment scores at 10%. In the following, we present results based on the union of upregulated (red) and downregulated (blue) proteins, as depicted in Figure 2 . GO term analysis was performed on up- and downregulated proteins separately, but these results were too similar to the union set, when combined, to warrant a separate discussion. The GO term enrichment analyses were conducted for both the control and clinical strains of P. aeruginosa ( Supplementary Figure S1 ).

Table 1A provides a list of the significantly enriched GO terms ( Supplementary Figure S1A ) with their respective statistics. Noteworthy terms include proteins associated with various cellular components (e.g., periplasmic space, cytosol, cytoplasm, and membrane), transporter activities (peptide transmembrane, efflux transmembrane, ABC-type sulfate and siderophore), DNA interactions, and specific processes (e.g., pyoverdine biosynthesis and iron ion transport). The clinical strain where only the fold-change threshold is applied, displayed in Table 1B ( Supplementary Figure S1D ), exhibits some of the same terms as in Table 1A , albeit with lesser significance overall. A unique feature here is the enrichment of proteins without associated GO term annotations, which, with an enrichment score of 1.52, is still significant due to the substantial number of proteins in the “NOGO” pseudo-term as GO annotations are missing.

To understand the proteomic landscape of imipenem resistance in P. aeruginosa, we have thus far delved into within-strain comparisons, specifically between the 1.0x MIC and control probe (no antibiotic) experiments. We have also cross-compared these results between the clinical and control strains. Now, we consider a more direct comparison across strains, by using the difference between the normalized, log2-transformed LFQ values of the control and clinical strains as a measure of fold-change (see Materials and Methods for details).

The volcano plots depicted in Figures 4A, C illustrate fold-change thresholds of 1.5 and 8 respectively (equivalent to approximately +/- 0.585 and exactly +/- 3 in log2 space respectively). The reason for these choices stems from Figure 2 , where we wanted a comparable example ( Figure 4A ) as well as a suitably sparse counter example, where only proteins with large differences in expression are included ( Figure 4C ). Other thresholds could be chosen, but we found a >8-fold change threshold to be a good compromise between sparsity and representativeness. The FDR corrected p-values displayed on the y-axis are derived from the control strain t-test. These are included for visual comparability with Figure 2 and are not considered in the GO term enrichment analysis displayed in Figures 4B, C .

As Figures 1B, C demonstrate, the clinical strain’s overall protein expression is not significantly affected by increasing MIC. Consequently, we anticipate that the normalized differential protein expression will mirror the results observed in the control strain. This expectation is confirmed in Figure 4B (and Supplementary Figure S1C ). A noteworthy change, however, is the enrichment of proteins associated with magnesium ion binding.

A GO term enrichment analysis was conducted using the 146 proteins with fold change of more than 8, denoted by blue and red dots in Figure 4C . The results of this analysis are presented in Figure 4C and further detailed in Table 1C . The proteins involved are provided in Supplementary Table S8 (in the “IDs” column).

Several of these terms are recognizable from protein subsets seen in e.g. Table 1A , but there are also differences from earlier protein subsets. These include transport terms like glycine betaine, amino-acid betaine, and carnitine transmembrane transport. There are specific enzyme-related terms such as sigma factor antagonist and cytochrome-c peroxidase, along with chaperone-mediated protein complex assembly. The list also includes cellular responses to hyperosmotic salinity and sulfur starvation, protein secretion by the type III secretion system, alginic acid biosynthesis, and structures like the extracellular region and host cell membrane. It is also important to highlight the n=33 proteins not associated with any known GO terms, representing a significant over-representation. These proteins, along with other specific proteins, will be the focus of our investigation in the subsequent section.

In Figure 4 , we visualized the differential protein expression comparing the difference in the clinical strain to that of the control strain, and in Figure 4C , we applied a threshold demanding this difference to be >8-fold. Here, we will delve into some of the proteins adhering to this threshold. To obtain a manageable number of proteins, we employed some criteria for selection: (i) low or negative correlation values (< 0.4) between clinical and control strain across the control probe (no antibiotic) plus 4 imipenem concentrations (0.1x, 0.25x, 0.5x and 1.0x MIC), (ii) striking differences in expression, and/or (iii) part of significantly enriched GO terms in Table 1C . This process led us to 26 proteins, 24 of which are shown in Figure 5 . The remaining two proteins are Tse4 (Q9I069), already displayed in Figure 3A , and PpiC2 (Q9HWK5) which is highly expressed in the clinical strain, and comparably non-expressed in the control across all concentrations. The remaining 24 proteins are split so that Figures 5A-O has an overall downregulation for the control for increasing imipenem concentrations, while Figures 5P-X show an overall up-regulation for the control. The clinical strain displays more stability across varying concentrations, which is not surprising given the results shown in Figures 1B, C . However, Figures 5S, W , alongside Figure 3A , displaying the proteins Wzm (Q9HTB8), PA0101 (Q9I731) and Tse4 (Q9I069) respectively, also show clear down-regulation in the clinical strain for increasing concentrations, in addition to the up-regulation in the control strain. As in Figure 3 , the standard deviations are shown as green (clinical) and orange (control) shading, where linear interpolation has been used between the measurements. While several proteins display large standard deviations for the clinical strain, most control strain trends are associated with comparably small standard deviations. This indicates a clear biological effect of increasing concentrations of imipenem for the P. aeruginosa (ATCC 27853) control strain which can be related to these proteins and their associated functions and compartments. More information about the proteins in Figure 5 can be found in Supplementary Table S2 , including their UniProt description (UniProt Consortium, 2023) and associated GO terms (Ashburner et al., 2000). Some proteins, PA2384 (Q9I195), PA4129 (Q9HWQ2), PA1469 (Q9HVV0), PA1034 (Q9I4U0), PA1404 (Q9I3U0) and PA0101 (Q9I731), in Figures 5C, D, I, M, P, W respectively, are uncharacterized and have thus no associated GO terms. An additional observation to remark in Figure 5 is the high relative expression of certain proteins in the clinical strain at 0 MIC (control probe), as well as the quasi-absence of other proteins, both groups largely displaying also the aforementioned stability over the different antibiotic concentrations. Both groups, constitutively high expressed or absent proteins in the clinical strain, belong largely to transporter protein and outer membrane and cell wall GO terms. These results will be further analyzed in the Discussion section, where we will shed light on their potential association with imipenem and antimicrobial activity.